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L  INTRODUCTION 

The  mativation  for  this  study  grew  from  interest  in  smart  materials  and  their 

appUcation  to  heticopter  dynamic  problems.  Smart  materials  can  be  used  for  aeroelastic 

tailoring  of  rotor  system  dynamic  components  and  also  as  a  mechanism  for  implementation 

of  active  rotor  control  without  adding  unwanted  mechanical  compl^^  Active  control 

of  helicopter  rotor  blades  has  several  potential  benefits.  Among  these  benefits  are 
reductions  in  vibration  and  acoustic  signature  as  well  as  elimination  of  mechanical  and 
aeromechanical  instabilities  such  as  air  and  ground  resonance.  The  idea  of  active  control 
of  heUcopter  rotors  has  recdved  much  attention  in  recent  years.  One  of  the  driving  forces 
behind  this  interest  is  the  increased  popularity  of  bearingless  rotors,  which  offer  the 
benefits  of  simplistic  design  and  maintainability,  but  often  give  rise  to  air  and  ground 

mechanical  stability  problems. 

It  was  immediately  apparent  that  in  order  to  conduct  an  in-depth  mvestigation  mto 
ftie  appfication  of  smart  material  technology  to  aeroelasticity  and  heficopter  active  control, 
the  ability  to  accurately  model  coupled  rotor-fiiselage  dynanncs  was  needed.  Several 
software  packages  such  as  CAMRAD®  by  Johnson  Aeronautics,  FUGHTLAB®  by 
Advanced  Rotorcraft  Technologies  Corporation,  and  UMARC®,  developed  at  the 
University  of  Maryland,  offer  this  capability.  All  of  these  codes,  while  very  capable,  are 
quite  complex  and  require  considerable  experience  in  order  fi)r  a  user  to  become  proficient 
in  applying  them.  As  a  result,  this  study  was  initiated  in  order  to  develop  a  computational 
tool  for  modeling  and  analysis  of  coupled  rotor  fiiselage  dynamics  by  utilizing  readily 
available  and  generally  applicable  technical  and  mathematical  software.  The  scope  of  this 
study  involves  utilizing  the  symbofic  manipulation  software,  MAPLE  ®  by  Waterioo 


1 


Software,  and  the  dynamic  emulation  software,  SIMULINK®  by  The  Mathworks,  to 
model  and  amulate  the  unstable  mechanical  phenomenon  of  helicopter  ground  resonance. 
The  resulting  models  will  be  used  to  do  parametric  studies  of  ground  resonance,  including 
a  look  at  the  effect  of  active  rotor  control  using  fuselage  state  feedback.  In  order  to  better 
quantify  the  effects  of  parameter  variations  and  active  control  schemes  from  amulation 
data,  the  moving  block  technique  was  utilized  for  time  history  analysis  to  quantify  modal 
damping  characteristics. 

Ground  resonance  is  a  potentially  destructive  mechanical  instability  that  can  occur 
in  helicopters  with  fully  articulated,  bearingless,  or  hingeless  main  rotor  designs.  The 
phenomenon  of  ground  resonance  is  the  result  of  a  coupling  between  fuselage  motion  on 
its  landing  gear  and  rotor  blade  lead-lag  and  flap  motion.  The  equations  of  motion 
describing  the  coupled  rotor-fuselage  system  are  nonlinear  and  generally  quite  complex 
even  for  sinq)lified  models.  Procedures  and  techniques  for  dealing  with  the  various 
complexities  of  ground  resonance  and  other  mechamcal  and  aeromechamcal  phenomena 
that  are  characteristic  of  helicopters  have  been  extenavely  investigated  over  the  past  few 
decades  and  an  abundance  of  literature  is  available  on  the  sulqect.  The  follovdng  section 
presents  a  brief  overview  of  various  approaches  and  techniques  used  in  modeling  coupled 
rotor-fiiselage  systems  for  the  purpose  of  studying  the  ground  resonance  problem  and 
active  rotor  control. 
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n.  BACKGROUND 


Hie  following  paragraphs  discuss  the  phenomenon  of  ground  resonance  and  some 
of  the  derivation  and  modeling  techniques  utilized  through  the  years  to  manage  the 
fonnulation  and  analysis  of  the  complex  equations  of  motion. 

Ground  resonance  has  been  an  observed  happening  in  lotorcraft  since  the  first 

autogyros  were  flown  early  in  the  20*^  century.  It  can  occur  when  any  rotor  system  is 

placed  on  a  flexible  support.  Essentially,  a  perturbation  of  a  rotor  blade  causes  the  rotor 
center  of  gravity  to  shift  creating  an  inertial  load  on  the  fuselage.  The  fuselage  is  flexibly 
connected  to  the  ground  via  its  landing  gear  and  will  start  to  oscillate  in  response  to  this 
inertial  load.  For  a  certain  range  of  rotor  rotational  speeds,  the  fuselage  oscillations  will 
cause  the  amplitude  of  the  blade  osdUations  about  their  respective  hinges  to  increase,  and 
this  further  increases  the  inertial  forces  on  the  fiiselage.  If  left  to  its  own  accord,  this 
coupling  of  fusdage  and  rotor  blade  motion  will  increase  in  amplitude  untfl  some  nonlinear 
restoring  force  brings  the  system  into  a  limit  cycle  or  until  some  part  of  the  rotorcraft  feils 

[Ref  1]. 

No  discussion  of  helicopter  ground  resonance  would  be  complete  without 
considering  the  classic  wori.  of  Coleman  and  Fdngold  [Ref  2]  completed  in  the  1940’s. 
Coleman  and  Feingold  successfully  analyzed  the  coupling  between  a  rotor  and  fuselage 
and  identified  the  ground  resonance  instability  as  a  purely  self  excited,  elasto-mechanical 
phenomenon.  Their  study  was  based  on  a  simplified  three-bladed  rotor  model  which  is  the 
basis  for  one  of  the  models  considered  in  this  study.  The  model  aUows  for  hub 
translational  d^ees  of  freedom  in  one  plane  and  rotor  blade  lead-lag  degrees  of  freedom 
in  the  same  plane  (see  Figure  3 . 1).  In  deriving  the  equations  of  motion,  Coleman  makes 


use  of  “Coleman”  coordinates  [Ref  2]  and  complex  variables  to  reduce  the  number  of 
equations  that  describe  the  system  to  two  complex  (four  real).  From  the  equations  of 
motion,  the  characteristic  equation  is  derived  by  assuming  a  solution  that  has  the  rotor 
center  of  mass  and  fuselage  center  of  mass  moving  in  an  elliptic  vdurhng  motion.  The 
roots  of  the  characteristic  eqi^on  are  the  characteristic  whirling  speeds  of  the  rotor ,  and 
the  nature  of  these  roots  indicate  the  nature  of  the  system  stability,  i.e.,  whether  or  not  the 
rotor  rotational  speed  resides  in  the  self  exdted  region.  Coleman  and  Fdngold  uWmatdy 
reduce  the  results  of  their  study  to  a  series  of  gr^hs  wduch  can  be  qiplied  to  a  wide  range 
of  rotor  configurations. 

Coleman  and  Feingold’s  work  became  the  basis  for  the  evolution  of  theory  and 
design  techniques  used  fijr  dealing  with  groimd  resonance.  Although  this  classic  theory 
oflfers  much  insight  and  understanding  into  the  phenomenon,  especially  for  conventional 
articulated  rotor  systems,  the  increasing  popularity  of  hingeless  and  bearingless  rotor 
designs  in  modem  helicopters  and  the  growing  desire  to  eliminate  the  need  for  mechamcal 
danqiers  requires  more  sophisticated  analytical  techniques. 

As  computational  power  improved  with  the  evolution  of  dig^  compirters,  more 
general  techniques  for  analyang  rotor  system  stability  came  into  being.  Peters  and 
Hohenemser  [Ref  3]  qiply  Floquet  analysis  to  the  problem  of  lifting  rotor  stability. 
Roquet  analysis  is  a  method  which  can  be  used  to  determiae  the  stability  of  solutions  to 
systems  of  linear  ordinary  differential  equations  with  periodic  coefficients.  The  Roquet 
transition  matrix  vdiich  relates  the  system  state  variables  at  the  beginning  and  end  of  a 
rotational  period  is  computed  by  numerical  time  wise  integration.  The  dgenvalues  of  the 
transition  matrix  are  a  measure  of  system  stability.  Hammond  [Ref.  4]  applies  Roquet 
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analysis  to  the  prediction  of  mechanical  instabilities,  specifically  examining  the  case  of 

lead-lag  damping.  The  unbalanced  problem  requires  solution  of  the  equations 
of  motion  with  the  blade  equations  expressed  in  the  rotating  reference  frame  because  a 
transformation  to  the  fixed  system  is  no  longer  pos^le  for  a  ground  resonance  anatyas  as 
was  pos^le  for  the  isometric  case.  As  a  result,  you  are  left  with  a  system  of  equations 
with  periodic  coefficients  which  can  be  handled  by  the  Floquet  method. 

Hingeless  and  beaiingjess  rotor  configurations  often  fece  the  additional  difficulty 
of  air  resonance.  Aerodynamics  may  play  more  than  a  passive  roll  in  the  ground 
resonance  regime  in  hingeless  ^sterns  in  contrast  to  articulated  ^sterns  wha"e 
aerodynamics  have  little  effect.  As  a  result,  more  complex  models  are  required  to 
accurately  represent  the  physics  of  the  heUcopter  aeromechanical  stability  problem. 

Modds  must  include  blade  fiap  and  torsional  degrees  of  freedom  as  well  as  lead-lag 
degrees  of  freedom.  Fuselage  models  also  should  include  pitch  and  roll  as  well  as 
translational  degrees  of  freedom.  Aerodynamic  models  can  range  from  quasi-steady  strip 
theory  to  unsteady  aerodynamic  theories  which  include  elaborate  wake  models  or  dynamic 
inflow  models.  Onniston  [Ref.  5]  utilizes  a  rigid  blade  and  rigid  fiiselage  model  with  flap- 
lag  and  pitch-roll  d^ees  of  freedom  to  conduct  parametric  investigations  based  on  an 
eigenvalue  analysis.  As  is  typical,  the  equations  of  motion  were  derived  by  a  Newtonian 
approach  and  the  resulting  system  of  nonlinear  differential  equations  are  linearized  for 
small  perturbations.  The  model  includes  linear  rotor  blade  and  landing  gear  springs, 
viscous  and  quaa-steady  aerodynamics.  Freidmann  and  Venkatesan  [Ref  6] 

and  Freidmann  and  Warmbrodt  [Ref  7]  derive  the  complete  set  of  governing  equations  of 
a  helicopter  rotor  coupled  to  a  rigid  body  fuselage.  The  equations  account  for  rotor  blade 


elastic  deformations  and  include  quasi-steady  aerodynamics  or  modified  Theordorsen 
unsteady  aerodynamic  theory.  In  deriving  the  full  equations  of  motion,  Freidmann  et  al., 
stress  the  importance  of  ^)plying  an  ordering  scheme  to  the  process  in  order  to  handle  the 
complexity  of  the  equations  and  enormous  number  of  terms  generated  by  their  expansion. 
The  equations,  as  presented  by  Freidmann  et  al.  [Ref  6  and  7],  are  in  a  form  which  makes 
them  generally  ^plicable  to  a  wide  range  of  rotorcrafl  problems. 

Another  interest  in  the  study  of  helicopter  ground  resonance  is  the  effect  that 
nonlinear  elastic  and  damping  forces  have  on  stability.  Tongue  [Ref  1],  Tongue  and 
Flowers  [Ref  8  and  9],  Tongue  and  Jankowski  [Ref  10],  and  Tang  and  Dowell  [Ref  11], 
use  variations  of  the  nonlinear  technique  of  harmonic  balance  uang  describing  functions  to 
r^esent  nonlinear  danqnng.  The  tedmique  is  useful  fi5r  investigating  limit  cycle 
behavior  of  strongly  nonlinear  systems  and  its  impact  on  system  stability. 

Active  control  of  rotor  systems  and  hs  qjplication  to  stabilizing  ground  and  air 
resonance  has  been  investigated  by  Straub  [Ref.  12]  and  Straub  and  Wambrodt  [Ref.  13]. 
In  both  of  these  studies  the  nonlinear  periodic  equations  of  motion,  derived  with  a  method 
■gimilar  to  that  of  Freidmann  and  Venkatesan  [Ref  6],  are  linearized  and  incorporated  into 
a  state  ^ace  model  in  which  active  control  inputs  are  iiq)ut  to  the  rotor  blades  fi'om  the 
fixed  coordinate  system  via  a  swaslqjlate.  The  state  space  model  is  then  used  to  study  the 
influence  that  state  feedback  gain  and  phase  have  on  system  damping. 

Helicopter  aeromechanical  instabilities  can  be  analyzed  by  methods  ranging  fi'om 
Coleman’s  classic  analysis  to  direct  time  integration  of  the  equations  of  motion.  As 
engineers  strive  to  develop  rotor  systems  free  of  ground  and  air  resonance  which  do  not 
require  the  addition  of  maintenance  iirtensive  mechanical  damping  systems,  more 
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elaborate  models  will  be  needed  to  accurately  cqmire  all  phydcal  aspects  of  the  problem. 
To  achieve  the  truly  daiiq)eriess  rotor  Ormiston  [Ref.  14]  addresses  three  differ^ 
approaches  which  may  be  feasible,  1)  incorporating  high  danqnng  material  into  the  blade 
or  fle?d)eam  structure,  2)  automatic  feedback  control,  and  3)  developm^  of  aeroelastic 
couplings  to  provide  inherent  stability.  These  three  approaches  have  provided  the  impetus 
behind  the  work  which  follows.  The  goals  of  the  following  study  were  to  develop  a 
mnHftling  technique  iitiliTing  ^mbolic  processing  to  manage  the  conq)lexity  of  deriving 
and  coding  the  coupled  rotor-fiiselage  equations  of  motion,  incorporate  the  resulting 
model  into  a  dynamic  simulation  environment,  and  have  a  final  product  whidi  provides  a 
useful  tool  for  conducting  parametric  studies  of  helicopter  am)mechanical  behavior. 
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ffl.  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  a  coupled  rotor-foselage  ^stem  were  derived  uring 
the  T  ^grangian  approadi.  This  study  was  conconed  with  two  models,  one  bdng 
characterized  as  simple  and  the  other  as  complex.  The  Lagrangian  ^proach  was  selected 
as  opposed  to  the  Newtonian  because  it  is  easily  implemented  with  the  aid  of  a  symbolic 
processor.  For  a  historical  note,  it  is  interesting  that  Lagrange  himself  recognized  the 
suitability  of  his  methods  to  routine  procesang.  He  states  in  YasM^dumique  Anafytique, 

The  methods  which  I  present  here  do  not  require  dther 
constructions  or  reasonings  of  geometncal  or  medi^cal  nature,  but  only 
algebraic  operations  proceeding  after  a  regular  and  uniform  plan.  Those 
who  love  the  Analysis,  will  see  with  pleasure  Mechanics  made  a  branch  of 
it  and  win  be  grateful  to  me  for  having  thus  extended  its  domain. 

The  equations  for  both  models  considered  were  formulated  with  Lagranges  method  in 

their  full  nonlinear  forms,  i.e.,  no  ordering  scheme,  small  angle  assumptions,  or 

^itipariTafinn  techniques  were  q^lied  during  derivation  and  subsequait  coding. 

A.  SIMPLIFIED  MODEL  DESCMPTION 

The  sinq)lified  model  is  based  on  that  used  by  Ck>leman  [Ref  2] ,  and  is  riiown 

schematically  in  Figure  3.1.  A  three  bladed  model  will  be  the  only  case  considered  in  this 

report,  but  all  mathematical  modeling  methods  used  in  this  study  can  be  easily  generalized 

to  any  number  of  blades.  Elastic  forces  generated  by  rotor  blade  and  flexbeam  motion 

were  modeled  as  a  linear  torsional  spring  located  at  the  effective  hinge  position  of  the 

blade.  The  gear  stiffeess  was  also  modeled  with  linear  springs.  For  the  baric 

simplified  model,  landing  gear  and  lead-lag  damping  was  modeled  with  linear  dashpot 

type  dampers.  Addition  of  nonlinear  mechanical  effects  sudi  as  hardening  and  softening 

springs,  hydraulic  damping,  and  lead-lag  stops,  wfll  be  discussed  in  a  later  section. 
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F^re  3. 1  Scheoutic  of  Simplified  Rotor-F^iseb^e  System 

This  model  allows  for  the  following  d^ees  of  freedom: 

Ui  =  Fuselage  translation  in  l-direction  (x-direction). 

U2  =  Fuselage  translation  in  2-direction  (y-direction). 

Ck=  Lead-lag  angular  displacement  of  k***  rotor  blade. 

A.  COMPLEX  MODEL  DESCRIPTION 

The  complex  model  is  based  on  that  used  by  Straub  [Ref  12].  This  model  assumes 
rigid  blades  and  frselage.  The  blade  and  fle7d}eam  elastic  forces  are  modeled,  as  in  the 
simplified  model,  as  equivalent  torsional  springs  located  at  effective  hinge  positions  ofi&et 
from  the  rotor  hub  (flap  and  lead-lag  hinges  are  assumed  to  be  coincident). 

This  model  allows  for  the  following  d^ees  of  freedom: 
ui  =  Fuselage  translation  in  1 -direction  (x-direction). 


U2=  Fuselage  translation  in  2-direction  O^-direction). 
ri  =  Fuselage  rotation  about  1-axis  (roll). 
r2= Fuselage  rotation  about  2-axis  (pitch). 

Lead-lag  angular  displacement  of  k***  rotor  blade. 

^ = Flap  angular  displacement  of  k*  rotor  blade. 

C.  COORDINATE  SYSTEMS  AND  TRANSFORMATIONS 

In  devdoping  the  equations  of  motion  for  the  two  coupled  rotor-fuselage  models 
five  coordinate  systems  were  utilized  with  transformations  between  the  various  systems 
based  on  Euler  angle  rotations.  The  five  coordinate  systems  are  (1)  inertial,  fixed  relative 
to  the  Earth,  (2)  fuselage,  fixed  to  center  of  gravity  of  fuselage,  (3)  hub,  parallel  to 
fiiselage  system  but  ofifeet  a  distance  h  in  the  positive  z  (or  3)  direction,  (4)  undefonned 
blade,  fixed  to  the  effective  hinge  position  on  the  k*  blade,  (5)  deformed  blade,  fixed  to 
the  e&ctive  hinge  position  on  the  k’^  blade,  but  with  the  x-axis  coincident  with  the  blade 
‘elastic’  axis.  Table  3.1  summarizes  the  notation  used  for  the  various  coordinate  ^ems. 


Table  3.1  Coordinate  Syston  Notation 


Coordinate  System 

Representative  Notation  for  a  Vector  in  tliis 
System 

Inertial 

y 

T 

Fuselage 

[y  y  z] 

T 

Hub 

X  y  Z 

Undefonned  blade 

f  /V  /V  ^ 

\x  y  z 

r 
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Coordinate  System 

Representative  Notation  for  a  Vector  in  this 

System 

Deformed  blade 

X  y 

T 

The  foDowing  generic  transformations  are  defined: 


1  0  0 
0  cos{a)  sin(a) 

0  -  sin(a)  cos(a) 


TSa)  = 


cos(a)  0  -  sin(a) 
0  1  0 

sin(a)  0  cos(or) 


TM  = 


cos(a) 
-  sinia) 


0 


sin(a)  0 
cos(a)  0 
0  1 


(3.1) 


(3.2) 


(3.3) 


Where  in  general,  Ti,  T2,  and  T3  are  rotations  about  the  1, 2,  and  3  axes  respectively 
This  notation  can  be  directly  utilized  with  the  ^bolic  processor  and  will  be  used  in  the 
following  section  where  the  energy  expressions  necessary  for  the  Lagran^an  derivation 
are  defined.  The  order  of  fuselage  rotations  when  using  these  Euler  angle  transformations 
win  be  pitch  -  roU,  and  the  order  of  rotor  blade  angular  displacements  wiU  be  flap-lag. 

The  Mowing  relations  summarize  the  coordinate  transformations  used  for  the  simplified 
and  complex  mod^; 

Simplified  model: 

1 .  Hub  to  Inertial: 
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;  Systems  are  parallel 


(3.4) 


X 

X 

’l  0  o' 

y 

=  I 

y 

;  /  = 

0  1  0 

z 

z 

O 

o 
_ 1 

2.  Blade  undefonned  to 


3 .  Blade  deformed  to  Blade  undefonned: 


Complex  model: 

1.  Fusdage  to  Inertial: 


2.  Hub  to  Fuselage: 


X 

X 

*1 

0 

o' 

y 

=  / 

y 

;  1  = 

0 

1 

0 

;  Systems  are  parallel. 

z 

z 

0 

0 

1 

(3.5) 


(3.6) 


(3.7) 


(3.8) 


3.  Blade  undeformed  to  Hub: 


13 


(3.9) 


4.  Blade  deformed  to  blade  xjndefbrmed; 


[j;te)3;(A)r 


(3.10) 


D.  DERIVATION  UTILIZING  SYMBOLIC  PROCESSOR 

This  section  summarizes  the  development  of  the  energy  expressions  necessary  for 
the  Lagrangian  derivation.  Here,  the  equations  for  the  complex  model  are  developed  to 
illustrate  how  the  symbolic  processor  was  utilized. 

The  Langrange  equation  can  be  expressed  as  follows: 


d  ( df\  ST  SU  SD 
dt\SqJ  Sq^^  Sqt  Sq^ 


(3.11) 


Where,  T  is  the  kinetic  energy,  U  is  the  potential  energy,  D  is  the  dissipation  fonction,  Fj  is 


a  generalized  force,  and  q;  is  a  generalized  displacement.  The  generalized  force  term,  Fi, 
will  describe  the  aerodynamic  forces  on  the  individual  rotor  blades  and  will  be  discussed  in 
a  later  section,  as  a  result,  this  derivation  develops  only  the  system  of  homogeneous 
equations.  The  various  energy  terms  can  be  broken  down  into  two  categories,  terms  due 
to  blade  motion  and  terms  due  to  foselage  motion,  to  ^e  the  following  equations; 

7-=7>+t(j;).  P.12) 

*=1 
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(3  13) 


ife=l 

/)=z)p+i;(c/,),  (3.14) 

i=:l 

Where  the  subscripts  F  and  B  indicate  fuselage  and  rotor  blade  re^>ectively. 

The  kinetic  energy  of  the  k*  rotor  blade  is  given  by  the  follo^ving  expresdon: 

(‘^BX=jj'"pp)dr  (3.15) 

Where  p  is  the  position  of  a  point  on  the  dastic  axis  of  the  k*  rotor  blade  with  respect  to 
the  inertial  coordinate  system  at  any  instant  in  time,  and  m '  is  the  mass  distribution  per 
unit  length  of  the  blade  (for  this  study  mass  distribution  per  unit  length  is  assumed  to  be 
uniform).  The  porition  of  a  point  on  the  elastic  axis  of  a  rotor  blade,  p ,  is  expressed  as 
the  sum  of  relative  positions  with  respect  to  the  various  coordinate  systems  transformed  to 
the  inertial  system.  Thus, 

P  —  {pFj)j  ,  (316) 

where,  for  example,  the  term  {pbu_h}^  is  the  portion  of  the  origin  of  the  undeformed 

blade  coordinate  system  with  respect  to  the  hub  coordinate  system  transformed  into  the 
inertial  coordinate  ^stem.  The  individual  terms  of  equation  (3.16),  referring  to  equations 
(3.1)  through  (3 .3)  are  defined  as  follows; 
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II 

1 

(3.17) 

(pf.f),  =\T,(r,h{n)J  Ph.f 

(3.18) 

(3.19) 

(ftM.*,),  =  =  0  Origins  are  camcident 

(3.20) 

(pr.ee),  =  h.ee 

(3.21) 

Where, 


«1 

’o' 

e 

r 

li 

“2 

0 

»  Ph_f  ~ 

0 

li 

0 

0 

>  Pp_Bd  ~ 

0 

0 

Equations  (3  .22)  are  substituted  into  the  equations  of  (3.17)  through  (3.21)  and 
the  matrix  multiplication  is  performed  with  the  results  substituted  into  equation  (3.16). 
This  gives  a  vector  expression  for  the  position  of  an  arbitrary  point  on  the  elastic  axis  of 
the  k*  rotor  blade  with  respect  to  the  inertial  coordinate  system  at  any  instant  in  time  in 
terms  of  the  system  degrees  of  freedom.  The  time  derivative  of  this  expression  gives  the 

velocity,  p ,  which  is  substituted  into  equation  (3.15)  to  give  the  kinetic  energy  for  k* 
rotor  blade.  All  of  the  calculus  and  algebra  was  accomplished  with  MAPLE®  (see 
Appendix  A  for  a  look  at  the  MAPLE®  worksheet). 

The  elastic  forces  generated  by  rotor  blade  motion  gjve  rise  to  a  potential  energy 
term  in  the  Lagrange  equation.  Since  a  rigid  blade  model  was  assumed,  the  potential 
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energy  was  modeled  uang  equivalent  toraonal  springs  to  restrain  the  rotor  blade,  with 
spring  constants  selected  to  approximate  elastic  forces  due  to  in  plane  and  out  of  plane 
bending  of  the  rotor  blade  (and  the  flexbeam  in  the  hingeless  case).  The  potential  energy 
of  the  k*  rotor  blade  is 

An  explanation  of  the  validity  of  using  an  equivalent  torsional  spring  system  to  model  the 
elastic  forces  of  a  deformed  rotor  blade  is  ^en  in  some  detail  by  Venkatesan  and 
Friedmann  (Kef  6]. 

System  ^  modeled  in  energy  form  by  use  of  a  disapation  fonction,  which 

for  the  k*  rotor  blade  of  the  complex  rotor  model  is 

For  the  fuselage,  the  kinetic  energy  in  terms  of  translational  and  rotational  d^ees 
of  freedom  is 


The  fiiselage  potential  energy  is 


(3.25) 

(3.26) 

is 

(3.27) 

(3.28) 

The  disapation  functions  for  the  fuselage  are 
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(D,L=|a!.r>|cM^  (3-30) 


The  resulting  inputs  for  the  fuselage  terms  in  equations  (3.12)  through  (3.14)  are 

■  (3.51) 


AD  of  the  energy  e5q)ressions  defined  above  were  entered  into  a  MAPLE® 
worksheet  programmed  to  apply  equation  (3.11)  and  generate  the  equations  of  motion 
corresponding  to  each  of  the  system’s  degrees  of  fi-eedom  (see  ^pendix  A).  An 
important  characteristic  of  MAPLE®  is  that  it  allows  for  distinction  between  dependent 
and  independent  variables  via  functional  notation,  e.g.,  to  indicate  a  variable  X’  is  a 
function  of  time  (t)  simply  write  it  as  ‘X(t)’.  It  is  also  important  to  note  that  when 
applying  Lagrange’s  equation  in  MAPLE®,  derivatives  are  only  understood  when  taken 
with  respect  to  independent  variables,  so  when  taking  derivatives  with  respect  to  the 
degrees  of  fi-eedom  and  the  time  rates  of  change  of  the  degrees  of  fi-eedom,  the  time 
functional  notation  which  represents  these  variables  must  be  converted  to  independent 
variable  notation.  For  example,  the  fl^  angle  degree  of  fi-eedom,  p  it) ,  and  its  time  rate 


of  change , 


spit) 

dt 


would  have  to  be  replaced  in  all  of  the  energy  expressions  by  the 


independent  variables  ,p  md  dp,  respectively,  in  order  for  terms  like  and 
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ilfwhere  ,,  =^«  and  «,  =^)  to  be  evaluated  pr<,periy  by  the  MAPLE* 

dqi 


to  be  evaluated 


-iflll  1 

symbolic  engine.  AddWonaUy,  for  the  time  derivative  term, 
psop^ty.dl  degree  ofifoedomeapreasedinindepm^ 

tacktolhnedependentnowion  The MAPI£* code udnchaeco«.p»estl«ata^^ 

aanipulalions  for  the  complex  model  is  contained  ta  An^ 

The  eiparions  of  motion  generated  by  the  MAPLE*  jnogram  for  the  simplified 

model  wereverifiedagainsttheequationsusedbynowersandTo^  FK— 

mn,  Tongue  also  utilized  a  Lagra.^  approach  and 

gaovving  equations  of  nKUionforanKuielsinnlar  to  the  one  d^^ 


M.y*C,y+r,M+^^yy  = 


"y  ^  y 

-muR 


sin(y.  +<.)-{Q+f.f  ein(r»  +4;)] 


m, 


R^L  +Ctik+  4*  +  m,a^^R  sm(4t)  - 

+4;)-)^cos(v^t  +4;)) 


(3.34) 


(3.35) 


(3.36) 


Flowers 


and  Tongue  included  an  additional  tenn  in  equations  (3.34)  and  (3.35)  mvotvmg 
am  pmduct  of  first  derivatives  and  their  absolute  values.  These  terms  r^  nonlinear 

damping.  Incfotionofoonlinear(hydiauHc)  damping  in  the  dissipati™ 

into  L^range’s  equation  will  be  discussed  in  a  later  section 
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Equations  (3.34)  through  (3.35)  were  compared  against  the  equations  of  motion 
generated  by  the  MAPLE®  program  for  a  three  bladed  simple  rotor  fuselage  model  and 
were  foimd  to  match  exactly  (except  for  the  nonlinear  damping  terms  Miich  were  not 
included  initially).  The  MAPLE®  equations  for  this  case  are  shown  in  Appendix  B. 
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IV.  BUILDmG  THE  SIMCIAHON  MODEL 


A.  S-FUNCnONS  AND  CODE  GENERATION 

Constructioii  of  the  simulation  model  from  the  equations  of  motion  was  based  on 
the  structure  of  the  SIMULINK*^  S-fimction.  The  S-fimction  defines  the  dynamics  of  a 
model.  It  can  be  written  in  C  or  Fortran  code  or  as  a  MATLAB®  m-file  (a  mathematical 
programming  language  with  similar  ^mtax  to  Fortran).  The  structure  of  the  S-fimction  is 
generic  so  as  to  allow  for  a  wide  range  of  fimctionality  programming  the  dynamics 

of  various  ^sterns.  SIMULINK*^  accesses  an  S-fimction  through  its  numerical  integration 
routines.  The  routines  make  calls  to  the  S-fimction  for  ^edfic  information,  the  type  of 
information  returned  is  dependent  on  the  value  of  a  flag  variable  sent  by  the  int^ration 
routine.  For  example, 

flag  =  0  S-fimction  returns  sizes  of  parameters  and  initial  conditions , 

flag  =  1  S-fimction  returns  state  derivatives  dx/dt, 
flag  =  3  S-fimction  returns  outputs. 

The  section  of  the  S-fimction  which  conqmtes  the  derivatives  at  each  time  step  is  the 
section  \\^ch  contains  the  equations  of  motion  [Ref  15]. 

As  a  resuh  of  the  Lagrangian  derivation ,  MAPLE®  generated  the  equations  of 
motion  in  the  following  form, 

F(^Xx,t)  =  Q  (4.1) 

where  x  is  a  vector  of  displacement  degrees  of  freedom  of  the  system.  Unfortunately 
this  form  is  not  very  useful  when  it  comes  to  programming  a  SIMULINK®  S-fimction,  so 
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MAPLE®  was  further  used  to  manipulate  the  equations  of  motion  into  the  following 
equivalent  form, 

[^(x,x,/)]  X  =  /(x,x,/)  (4.2) 

where  is  an  matrix  and  /  an  Nxl  vector,  with  N  =  number  of  degrees  of  freedom 
of  the  system.  This  is  possible  since  the  equations  are  quasi-linear  in  the  second  derivative 
(acceleration)  terms,  i.e.,  no  terms  exist  of  types  such  as  x%  or  sin(x) ,  etc.  This  form 
can  then  be  transformed  from  N  second  order  equations  to  2N  first  order  equations  as 
follows. 


x  =  w 


(4.3) 


These  equations  can  be  evaluated  at  each  time  step  in  a  numerical  simulation  to  give  the 
state  derivatives.  The  primary  job  to  be  accomplished  with  MAPLE®  was  to  generate  the 
expressions  for  the  elements  of  [A]  and  /  (\^ch  can  be  quite  lengthy)  from  the 
equations  of  motion.  After  this  was  accomplished,  the  MAPLE®  code  generation  routine 
was  used  to  automatically  generate  the  optimized  C  or  Fortran  code  that  could  be  placed 
directly  into  an  S-fimction  template  (See  Appendbc  C  for  an  example  of  the  MAPLE® 
code  generation  resuhs). 

For  each  model,  MAPLE®  was  used  to  generate  the  Fortran  code  representing  the 
equations  of  motion  This  code  was  used  to  create  an  S-fimction  in  MATLAB®  m-file 
format.  Because  Fortran  and  m-file  syntax  are  so  sinular,  only  minor  editing  was  required. 
A  copy  of  the  S-fimction  program  for  this  case  is  contained  in  Appendix  D. 
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B.  NUMERICAL  INTEGRATION  ROUTINE  (ODE  SOLVER) 

SIMULINK®  provides  several  numerical  ordinary  diflFerential  equation  (ode) 

solvers  (nimierical  integration  algorithim  The  algorithnis  utilized  in  this  study  are  from 

the  Runge-Kutta  (rk)  femily  (ik45  and  rk23).  Runge-Kutta  algorithms  generally 
outperform  other  schemes  for  systems  of  nonlinear  ordinary  differential  equations  which 
are  not  too  stiff.  The  rk  algorithms  also  handle  discontinuities  well  [Ref.  7]. 

For  completeness  and  to  describe  how  the  int^ation  algorithms  intmct  with 
SIMULINK®  S-functions,  a  brief  description  of  the  rk  method  is  given.  The  first  order 
ode  (the  following  algorithm  can  be  easily  extended  to  systems  of  first  order  ode’s). 


can  be  integrated  between  and  t„+i  to  give  the  following. 


(4.5) 


In  the  rk  method,  the  integral  eiqiression  on  the  right  hand  side  of  equation  (4.5)  is 
approximated  with  a  numerical  integration  scheme,  such  as  Simpson’s  1/3  rule,  resulting  m 


the  following  eiqiresaon  for  the  value  of  at  the  next  time  step. 


(4.6) 


The  parameter  h  is  the  ^e  of  the  time  step,  and 


and  7^1  are  initial  estimates  of  at 


the  half  and  full  step,  and  are  given  by  the  following  expressions  derived  from  Taylor 


expan^ons  about  y,, , 
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(4.7) 


y  t 

n+- 

2 


h 

2 


y  18^  1 

(4.8) 


The  variable  ^  is  a  weighting  parameter  that  can  be  selected  to  optimize  the  accuracy  of 
the  numerical  method  [Ref  16],  In  SIMULINK®,  h  is  the  values  of  the  function  /(yjj 
that  are  generated  by  the  S-function  when  the  integration  routine  calls  with  the  proper 
flag.  Thus,  SIMULINK®  provides  the  means  of  modeling  any  dynamic  system  for 
numerical  simulation  provided  that  the  equations  of  motion  are  expressed  as  a  system  of 
first  order  ode’s.  It  is  important  to  note  that  the  algorithm  described  above  is  only  one 
variant  of  the  ik  family  of  ode  solvers  and  that  the  actual  routines  utilized  by 
SIMULINK®  are  somewhat  more  sophisticated. 

Before  he^nning  a  simulation  with  SIMULINK®,  the  user  sets  several  parameters 
which  control  the  execution  of  integration  routine.  The  user  must  designate  a  maximum 
and  minimum  (time)  step  size,  emulation  duration  (start  and  stop  time),  and  a  tolerance 
which  establishes  the  maximum  relative  error.  If  the  algorithm  cannot  decrease  the 
maximum  relative  error  without  going  below  the  minimum  step  size,  a  warning  message  is 
displayed. 
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V.  MODELING  NONLINEAR  EFFECTS 
To  include  additional  nonlinear  terms  into  the  overall  model  the  equivalent  energy 
expression  representing  the  effect  is  simply  added  to  the  overall  energy  expression  in  the 
MAPLE®  worksheet  program.  For  example,  to  model  a  nonlinear  flexbeam,  a  DuflBng 
spring  term  can  be  added  to  the  equation  of  motion  of  the  k*  rotor  blade  of  the  form 


For  flie  Lagrangian  derivatioii,  wUdi  the  MAPLE®  program  perfi>niis,  the  equivalent 
potential  energy  term  is  given  by 

This  term  is  amply  added  to  the  expression  represrating  the  potential  energy  of  the  k* 
rotor  blade.  To  add  nonlinear  damping  to  the  rotor  blades  or  fuselage  a  term  of  the 
following  form  is  added  to  the  respective  dissipation  function 

(5-3) 


where  Vx  is  the  nonlinear  hydraulic  damping  co^dent  [Ref.  5]. 

The  effect  of  lead-lag  (flap)  stops  was  also  modeled  by  incorporating  a  Emulated 
jump  in  lead-lag  (flap)  stiffiiess  by  use  of  Heaviside  step  fimctions  in  the  potential  energy 
e5q)resdons.  For  lead*lag  stops  the  expression  is 

-^k,h(c+z)(c+zY  (5.4) 


where  the  Heaviside  step  function,  H{r),  is  defined  as 

fO,  T<0 


(5.5) 
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and  z  is  the  absolute  displacement  angle  at  which  a  rotor  blade  would  engage  the  lead-lag 
(flap)  stops,  and  K,  is  the  effective  increase  in  lead-lag  (flap)  stififiiess.  For  programming 
puiposes,  e7q)ression  (5.4)  can  be  written  in  terms  of  the  signum(  )  function  as  follows 


z)  +  A:^^5(gw«777(^+z)  -  z) 

+  ^K,z^sigmmi(^-2}-jK^C  signum(C-2)-jK,C^sigmm(C+z) 


(5.6) 


where 


sigmm(x)  =  -^  (5.7) 

Thus,  any  structural  or  damping  nonlinearity  can  be  incorporated  into  the  model  by 
adding  the  appropriate  energy  expression  into  the  MAPLE®  worksheet,  and  then 
executing  the  worksheet  to  generate  the  updated  code  for  incorporation  into  the 
SIMULINK®  S-function. 
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VL  SIMULATION  RESULTS 

This  section  displays  results  of  several  simulations  and  demonstrates  the  unique 
capabilities  and  flexibility  of  the  modeling  method  described  in  previous  sections.  Direct 
simulation  allows  analysis  of  any  number  of  different  configurations  or  scenarios,  such  as 
non-isotropic  hub,  one  damper  inoperative,  or  even  simulated  rotor  blade  damage. 
Though  the  time  history  plots  in  the  following  subsections  do  not  indicate  h,  SIMULINK 
offers  the  useful  capaMity  of  bdng  able  to  visualize  tiie  dynamics  of  a  model  as  th^r 
progress,  which  can  add  valuable  msigbt  into  the  phenomenon  being  studied. 


A.  SIMPLE  MODEL 

The  following  table  summarizes  the  parameters  that  can  be  set  interactivety  for  any 
simulation  for  the  ample  model.  The  table  gives  the  representative  nomenclature  used  to 
represent  the  parameter  in  the  MAPLE®  and  S-function  programs.  The  table  is 


representative  of  a  three  bladed  model. 


Table  6.1  Simple  Rotor  Model  Pn^ram  Nmneadatore 


Paramrter 

As  h  Appears  in  MAPLE®  and/or 
S-function  code 

Units 

Rotor  blade  mass 

mbfn,mb(2),mb(3) 

mass 

Fuselage  effective 
mass  in  x  and  y 
direction 

M(1),M(2) 

mass 

Distance  fi'om  hinge 
to  center  of  mass  of 
blade 

R 

length 

Rotor  speed 

Omega 

rad/sec 

el 

length 

z 

radians 

Azimuth  phase  angle 
of  rotor  blade 

Phi(l),Phi(2),Phi(3) 

radians 
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Parameter 

As  it  Appears  in  MAPLE*  and/or 
S-function  code 

Units 

Lead-lag  linear 
damping  coefficient 

Czeta(l),  Czeta(2),  Czeta(3) 

mssi 

Vzeta(l),  Vzeta(2),  Vzeta(3) 

Fuselage  linear 
damping  coefficient 
in  X  and  y  direction 

c(l),  c(2) 

force/ 

(length/sec) 

Fuselage  nonlinear 
damping  coeffident 
in  X  and  y  direction 

v(l),  v(2) 

force/ 

(length/sec)^ 

Lead-lag  linear 
spring  coeffident 

Ke(l),Ke(2),Ke(3) 

moment/rad 

Lead-lag  nonlinear 
spring  coeffident 
(Duffing  spring) 

Kd(l),  Kd(2X  Kd(3) 

moment/rad^ 

Lead-lag  stop  spring 
coeffident 

Ks(l),  Ks(2),  Ks(3) 

moment/rad 

Effective  fuselage 
stiffiiess  in  x  and  y 
direction 

K(1),K(2) 

force/length 

Fuselage  states 
initial  displacement 
conc^ons 

xXi,xYi 

length 

Fuselage  states 
initial  rate  conditions 

xrXi,  xrYi 

length/sec 

Blade  states  initial 
displacement 
conditions 

xli,  x2i,  x3i 

rad 

Blade  states  initial 
rate  conditions 


xrli,  xr2i,  xr3i 


rad/sec 


Table  6.2  shows  the  baac  simulation  case  for  the  parameters  in  table  6. 1 .  These 
base  values  will  serve  as  a  starting  point  for  each  simulation,  i.e.  \n^en  a  spedfic  parameter 
is  changed  for  a  amulation,  it  is  understood  that  all  other  parameters  will  be  set  to  the 
base  case  values. 


Table  6J,  Parameter  Settings  for  Basic  Simulation  Case 


mmmmm 

el 

z 

1  lOft 

170  radians/sec 

0.5  ft 

n/12  radians 

Phid) 

Phif2) 

Phi(3) 

0  radians 

2x/3  radians 

4%/3  radians 

c(l) 

c(2) _ 

HHKSIflHIi 

0  Ibs/fbs 

( 

3ibs/fps 

ibeib^shi 

Czeta(l) 

Czeta(2) 

Czeta(3) 

0  ft-Ibs/(radians/sec) 

0  ft-lbs/(radians/sec) 

0  ft-lbs/(radians/sec) 

Vzeta(l) 

Vzeta(2) 

Vzeta(3) 

0  ft-Ibs/(radians/sec)^2 

0  ft-lbs/(radian/sec)^2 

0  ft-lbs/(radian/sec)^2 

Ke(l) 

Ke(2) 

Ke(3) 

0  ft-lbs/radian 

0  ft-lbs/radian 

0  ft-lbs/radian 

Kd(l) 

Kd(2) 

Kd(3) 

0  ft-Ibs/radian'^S 

0  ft-lbs/radian^3 

0  ft-lbs/radian^3 

Ks(l) 

Ks(2) 

Ks(3) 

0  ft-Ibs/radian 

0  ft-lbs/radian 

0  ft-lbs/radian 

K(l) 

_ 5(2} _ 

113.000  Ibs/ft 

113,000  Ibs/ft 

xXi 

xTi 

xli 

x2i 

x3i 

Oft 

Oft  n 

Oradians 

Oradians 

Oradians 

xrXi 

xrYi 

xrli 

xi2i 

xr3i 

0.5  ft/sec 

Oft/sec 

0  radians/sec 

0  radians/sec 

0  radians/sec 

The  baac  case  is  intentionally  set  up  with  zero  damping  and  vnth  a  rotor  speed  set 
approximately  at  the  center  of  the  r^esang  lead-lag  mode  instability  r^on.  The  fira  set 
of  simulations  will  demonstrate  the  syaem  behavior  when  exdted  with  an  initial  fuselage 
velodty  as  indicated  in  Table  6.2.  Figure  6. 1  and  Figure  6.2  show  the  lead-lag  time 
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histories  and  the  fuselage  center  of  mass  trajectory  (displacements  are  in  feet)  for  the  basic 


case. 


iTigitra  6.1  Rotor  Lead-lag  Displacements  for  Basic  Parameter  Case  Settings,  Center  of  Self  Excited 

R^on. 

Unless  otherwise  specified,  for  plots  of  rotor  blade  motion,  i.e.,  lead-lag  and  flap , 
the  colors  red,  blue,  and  green  distinguish  the  individual  blades. 
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Figure  6,2  F^selj^e  Trajectory  for  Basic  Parameter  Settings,  Center  of  Seif  Excited  Region. 


As  expected,  Figures  6, 1  and  6.2  show  the  rapid  divCTgence  of  the  model  as  a 
result  of  being  in  the  center  of  the  self  excited  region  and  perturbed  with  an  initial 
fuselage  velocity  in  the  x-direction.  The  diverging  spiral  path  of  the  fuselage  center  of 
magfi  is  a  characteristic  result  of  the  regressing  lead-lag  mode  instability. 

Figure  6.3  and  6.4  show  the  corresponding  results  for  operation  below  the  self 
excited  region.  Figures  6.3  shows  a  beat  or  modulation  of  the  blade  response  but  no 
divergence.  The  beat  phenomenon  indicates  the  blade  lead-lag  motion  consists  of  two 
Hniwiinaiit  modes  closely  spaced  in  frequency.  The  fuselage  center  of  mass  trajectory 
shown  in  Figure  6.4  shows  an  elliptical  path  with  the  major  axis  of  the  ellipse  rotating 
about  the  zero  displacement  position.  Both  the  beat  phenomenon  and  the  precession  type 
motion  of  the  fixselage  seem  to  be  a  characteristic  behavior  of  a  system  operating  outtide 
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the  self  excited  region.  It  is  interesting  to  point  out  that  this  behavior  is  also  characteristic 


Figure  Rotor  Leai-lag  Timie  EQstories,  Rotor  Speed  Below  Self  Excited  Regioo 
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Trajectory  of  Fuselage  Center  of  Mass 


Figure  6.4  Fiiselage  Trajectory  for  Basic  Parameter  Settings,  Rotor  Speed  Below  Self  Excited 

Region 

Figures  6.5  and  6.6  show  the  results  of  a  simulation  where  rotor  speed  was  set 
above  the  self  excited  region.  Again,  the  fuselage  exhibits  an  elliptic  whirling  motion  with 
the  major  axis  of  the  ellipse  rotating  about  the  zero  displacement  position  while  the  blade 
lead  lag  motion  follows  a  beat  patt^. 
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YAxis  I  I  radians 


Figure  6.5  Rotor  Lead-lag  Time  ffistories,  Rotor  Speed  Above  Self  Excited  Region 


Figure  6.6  Fuselage  Center  of  Mass  Trajectory,  Rotor  Speed  Above  Self  Exdted  Region 
Figure  6.7  is  the  Coleman  stability  plot  [Ref.  2]  for  the  basic  configuration.  The 


red  lines  indicate  the  boundaries  of  the  self  excited  region  and  the  blue  line  marks  the 


center  of  the  self  excited  region.  The  X’s  indicate  the  operating  points  for  the  three  cases 
shown  in  Figures  6.1  through  6.6. 


Coleman  Stability  Chart  (Isotropic  hub) 


Mgare  6.7  Coleman  Stability  Plot  for  Basic  Case 


At  this  point,  a  comparison  was  made  between  the  simulation  model  and  a  time 
history  solution  of  Coleman’s  and  Feingold’s  equations.  Bramwell  [Ref  17]  derives 
Coleman’s  and  Feingold’s  equation  in  a  form  equivalent  to  that  of  the  simulation  model 
witih  the  blade  displacements  expressed  in  the  rotating  coordinate  system  and  the  fiiselage 
displacements  ejqjressed  in  the  fixed  coordinate  system.  These  equations  were  solved  in 
the  fixed  coordinate  system  using  an  dgenvalue  analysis  and  the  solutions  transformed 
back  to  rotating  coordinate  Systran.  A  comparison  was  then  made  with  the  lead-lag 
displacement  time  history  of  the  simulation  model.  Figure  6.8  shows  the  result  of  the 
comparison  using  the  parameters  of  the  bade  configuration  with  a  moderate  amount  of 
damping  added  to  rotor  blades  and  fuselage.  Figure  6.8  shows  excellent  agreement 
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between  the  two  solutions  with  a  significant  departure  between  the  two  occurring  only 
when  displacements  get  very  large.  Thus,  for  the  limiting  case  of  an  isotropic  hub  with 
linear  spring  stiffiiess  and  damping,  the  above  comparison  offers  some  amount  of 
verification  as  to  the  accuracy  of  the  simulation  model. 


Moving  on  from  the  basic  results  and  model  verification,  some  of  the  more 
int^esting  cases  that  were  simulated  will  now  be  discussed.  Figure  6.9  shows  a 
comparison  between  a  case  where  aU  blade  lead-lag  dampers  are  operating  and  a  case 
where  one  damper  is  inoperative.  The  first  plot  of  Figure  6.9  shows  a  rotor  with  all  blade 
dampers  operating,  in  the  second  plot,  the  blue  damper  is  disabled  by  reducing  the 
Hampitig  coefficient  by  two-thirds.  As  is  evident  from  the  plot,  the  very  sli^tly  unstable 
case  with  fiill  damper  operation  is  made  highly  unstable  by  fiiiling  one  damper. 
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Figure  6,9  Ooe  Lead-Lag  Damper  Dioperative 
Figure  6. 10  shows  the  results  of  simulating  damage  to  a  rotor  blade  by  reducing 

the  mass  of  the  blue  blade  by  20%.  The  undamaged  blades  are  forced  to  oscillate  around 

a  non-zero  displacement  position  in  order  to  compensate  for  the  damaged  blade,  but  the 

attiplitiniftH  of  all  the  blade  oscillations  appear  to  be  constrained. 
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Lea(Magdisp(rad)  ^  Lead-lag  disp  (rad) 


Lead-lag  t'me  histoi1es:sinnulated  one  blade  damaged 


0.1 


-0.1 

0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9 

lime  (sec) 

Figure  6.10  Simulated  One  Rotor  Blade  Damaged 
Figure  6.1 1  shows  the  effect  of  enabling  lead-lag  stops  in  the  model.  The  figure 
compares  the  time  history  of  a  blade  with  no  stops  with  that  of  a  blade  with  stops 
simulated  at  ±15  degrees.  Figure  6.12  shows  the  corresponding  fuselage  displacements. 
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y  displacement  |  I  lead-lag  displacement  (radians) 


Figure  6.11  Effect  of  Lead-Lag  Stops 


Figure  6.12  Fuselage  Displacements  •with  and  without  Simulated  Lead-lag  Stops 
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The  oiiyective  of  the  next  set  of  shmilations  was  to  examine  the  effect  of  a 
nonlinear  £I«d)eam  incorporated  Into  a  beaiin^ess  rotor  dedgn.  The  nonlinear  behavior  of 
the  fle^eam  was  assumed  to  be  that  of  a  DufBng  spring  where  the  restoring  moment  is 
given  by 

<*•*) 

is  the  linear  stiffeess  and  K4  the  nonlinear  stiffeess.  Simulations  were  conducted  for 
several  values  of  the  noirlinear  spring  constant  keeping  the  linear  coefhdent  constant  at 
22,000  ft-lbs/radian.  Results  are  shown  in  Figure  6.13.  The  primary  effect  of  increasing 
the  nonlinear  ^ring  constant  is  in  the  limiting  of  the  amplitude  of  the  lead-lag  response. 

As  is  depicted  in  the  Hgure  6. 13,  the  case  for  0  is  very  unstable  and  a  helicopter 

cau^  in  groimd  resonance  in  such  a  configuration  would  most  likdy  expoience 
catastrophic  feilure.  By  adding  the  hardening  (cutnc)  term,  the  unbovmded  growth  in 
amplitude  can  be  checked  as  is  apparent  fi'om  the  responses  for  the  cases  of  A!^=4£+5  and 
.^^8E+5.  As  the  anqrlitude  increases,  the  magnitude  of  the  nonlinear  term  becomes 
more  influential  and  effectively  changes  the  fi^equeiu^  of  osdllation,  shifting  it  outade  of 
the  unstable  r^on  and  allowing  the  osdllations  to  decay.  Once  the  amplitude  decays  to 
where  the  influence  of  the  nonlinear  term  becomes  snrall  the  cycle  repeats  itself  While  the 
limiting  an^litudes  for  the  noirlinear  cases  of  Figure  6. 13  are  still  large  for  lead-lag 
displacements  (on  the  order  of  30  to  40  degrees),  this  limiting  behavior  may  be  trough  to 
prevent  destruction  of  an  aircraft  if  groimd  resonance  were  to  be  excited.  In  f!i^ ,  whra 
lead-lag  displacements  are  small,  the  hardening  effect  of  a  nonlinear  fiexbeam  would  be 
negligible,  and  could  be  deagned  to  act  as  soft  inplane  in  order  minimize  hub  moments. 
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Figure  6.13  Effect  of  Hardening  Duffing  Flexbeam  on  Lead-lag  Response 


It  is  important  to  note  that  the  elastic  behaviors  of  tiie  flexbeams  modeled  by  the 
curves  in  the  upper  plot  of  Figure  6. 13  are  purely  hypothetical  and  were  selected 
arbitrarily  in  order  to  illustrate  the  effect  that  nonlinear  elastic  behavior  could  have  on 
rotor  system  response  and  stability. 


41 


B.  COMPLEX  MODEL 


The  complex  model  is  based  on  a  configuration  used  by  Straub  ^ef.  12],  The 
conq)uter  code  nomenclature  and  parameter  values  for  the  basic  case  used  to  conduct 
Emulations  for  this  study  is  contained  in  ^pendix  F.  .^)pendix  F  is  an  example  of  a  the 
MATLAB®  input  file  used  for  the  complex  model.  This  method  of  input,  as  opposed  to 
the  graphical  inteifoce  masking  feature  used  for  the  single  model  [Ref  15],  was  more 
convenient  in  the  case  of  the  complex  model  due  to  the  large  number  of  parameters. 

What  follows  are  examples  of  some  of  the  time  histories  generated  from  simulations 
conq)leted  with  the  complex  modd. 

Figures  6. 14  and  6. 15  show  the  fl^  and  lead-lag  response  of  the  rotor  to  a 
fuselage  roll  perturbation  (initial  angular  displacement  about  the  fuselage  x-axis).  For  this 
case  the  ground  resonance  was  not  excited  and  both  the  flap  and  lead-lag  motions  settle 
very  quickly.  Notice  that  the  lead-lag  di^lacement  settles  around  a  non-zero  steady  state 
poEdon.  This  is  due  to  the  aerodynamic  drag  on  the  rotor  blade,  the  modeling  of  wfaidi  is 
discussed  in  the  nect  section. 

The  next  set  of  figures  show  the  results  for  the  same  configuration  used  for 
Figures  6.14  and  6.15,  but  the  rotor  rotational  speed  has  been  changed  so  that  ground 
resonance  is  excited  by  the  roll  perturbation.  Figures  6. 16  and  6. 1 7  show  the  and 

lead-lag  response  for  this  case,  and  Figure  6. 18  shows  the  trajectory  of  fiiselage  pitch  and 
roll  displacements 
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Lead-^g  displacement  (radians) 


Figure  6.14  Flap  Response  to  Fuselage  Roll  Perturbation 


Figure  6.15  Lead-lag  Response  to  Fuselage  Roll  Perturbation 


43 


Lead-tag  displacement  (radians)  I  Flap  displacement  (radians) 


Figure  6.16  Flap  Response  with  Ground  Resonance  Excited 


Lead-Lag  Response  to  Fuselage  Roll  Perturbation  (Ground  Resonance) 


F^re  6.17  Lead-lag  Response  with  Ground  Resonance  Excited 
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Figure  6.18  Fuselage  Displacement  Trajectory  with  Ground  Resonance  Excited 


The  next  result,  shown  in  Figure  6. 19,  is  the  lead-lag  response  with  one  damper 
inoperative.  For  this  case,  the  system  has  enough  inh^ent  stability  to  settle  out  after  the 
initial  fiiselage  displacement.  The  undamped  blade  (blue  blade)  simply  settles  out  at  a 
higher  ampIttnHe^  but  this  anq)litude  is  small  enough  that  the  inertial  forces  that  arise  fi*om 
the  imbalance  have  a  small  effect. 
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Lead-lag  displacement  (radians) 


Figure  6.19  Lead*lag  Response  with  One  Damper  Inoperative 
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Vn.  INTRODUCING  AERODYNAMICS  TO  THE  MODEL 


Aerodynamic  forces  were  derived  using  quasi-steady  str^  theory.  Stall, 
con^essibOity,  reversed  flow,  and  wake  effects  are  ignored,  and  induced  flow  is  obtained 
ffom  momentum  theory.  The  results  of  the  derivation  are  the  aerodynamic  moments  about 
the  blade  hinge  in  the  flap  and  lead-lag  directions.  These  moments  are  ent^ed  as 
generalized  aerodynamic  forces  m  the  Lagtangian  derivation.  The  development  of  the 
aerodynamic  equations  was  adapted  Ifrom  the  approach  utilized  by  Kaza  and  Kvatemik 
[Rdf.  20]  for  a  fliap-lag  stability  analysis  on  rigid  articulated  blades.  The  main  difference 
between  thdr  approach  and  the  approach  outlined  in  the  following  paragr^hs  is  that 
fiiselage  motion  and  its  influence  on  blade  motion  is  contidered  here. 

The  relative  air  velodty  at  a  point  on  the  k*^  rotor  blade  due  to  forward  flight  at 
any  instant  in  time,  e^spressed  in  inertial  coordinates  is 

=  //  ORe.+Oe^  +XQRe,  (7.1) 

The  total  velodty  rdative  to  a  blade  element  is  obtained  if  the  above  expression  is  added 
to  the  n^ative  of  the  time  derivative  of  the  instantaneous  podtion  of  the  blade  element.  In 
inertial  coordinates  the  total  velodty  is  given  as 

(’-2) 

This  velodty  is  expressed  in  blade  deformed  coordinates  throu^  the  following 
transformation, 

= 3;fc)n(A)r,(n)2;(/-,)3;(r,X»L,),  (7.3) 

where  Ti,  T2,  and  Ts  are  given  by  equations  (3.1)  to  (3.3).  can  also  be  expressed 

as 
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UjCy  U  pi,  (7.4) 

Thus,  from  the  components  of  the  vector  ejqjression  (7.4),  the  radial,  tangential,  and 
perpendicular  components  of  the  total  velocity,  with  respect  to  blade  deformed 


coordinates,  are  given. 

The  lift  and  drag  acting  on  an  elemental  section  of  blade  are 

dL  =  —pacU'^adr 


dD  =  \-paclJ^ 


(C, 


<lo 


V  a 


dr 


where 


(7.5) 

(7.6) 


u  =  ^u/+u/ 


(7.7) 


The  angle  of  attack,  a,  is 


cr  =  ^-tan  ' 


0—  (ft 


(7.8) 


where  ^  is  the  section  pitch  angle.  The  lift  and  drag  are  then  transformed  to  give  the 
resultant  forces  along  the  y  and  z  axes  of  the  blade  deformed  coordinate  system,  thug 
giving 

dF.  =  -dL  sin(^  )  -dDoo^(f>  )  (7.9) 

dF-  =<2Lcos(^  )-ia[Dsin(^)  (7.10) 

To  obtain  the  generalized  aerodynamic  forces  from  the  above  expressions,  the  principle  of 
virtual  work  is  applied  for  a  fl^  -  lag  blade  displacement  sequence.  To  accomplish  this 
the  blade  differential  forces  given  by  (7.9)  and  (7.10)  are  transformed  to  the  blade 
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undefonned  coordinate  system  since  the  generalized  blade  displacements,  fiij(t)  and  ^t), 
are  e^qn-essed  in  this  fiame  of  reference.  The  transformation  is  as  follows 

[  0  ' 

(<<?)»,  =  dF,  (7.11) 

dF,\  \_dF,_ 

Applying  the  principle  of  virtual  woit,  the  generalized  aerodynamic  forces  on  the  k  rotor 
blade  are 

K). =f  <’  ‘2) 

K). = r  (“^L  •  = //<*>  (">3) 

Tsu^ere  is  the  position  vector  of  an  arbitrary  point  on  the  deformed  rotor  blade 

elastic  axis  with  respect  to  the  blade  undeformed  coordinate  system. 

To  simplify  the  incluaon  of  aerodynamics  into  the  model,  the  int^ral  e^qrressions 
for  the  generalized  aerodynamic  forces  appearing  in  equations  (7. 12)  and  (7. 13)  will  be 
evaluated  by  assuming  the  mean  value  of  the  forces,  dFp  and  dF^ ,  occur  at  the  r  =  0.7/? 

radial  position,  and  that  this  radial  position  also  corresponds  to  the  center  of  lift  and  drag 
on  a  rotor  blade.  With  these  Amplifying  assumptions,  the  resulting  generalized 
aerodynamic  forces  are 

The  inflow  ratio,  A ,  and  advance  ratio,  n  spearing  in  equation  (7.1)  are 
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(7.16) 


A  = 


(v  sin  or -v.) 


OR 


P  = 


V  cosa 
QR 


(7.17) 


The  induced  flow,  v,  ,  in  equation  (7.16)  is  calculated  by  equating  the  integrated  thrust  to 
the  thrust  from  momentum  theory,  leading  to  the  result. 


Cj.SiR 


(7.18) 


The  thrust  coeflBcient,  Ct  ,  is  determined  by  adding  the  average  total  lift  generated  by  each 
rotor  blade  over  one  rotor  revolution  and  dividing  the  quantity,  p  jt  Q^R* . 

For  this  study,  only  a  rotor-fuselage  ^stem  in  the  groimd  resonance  regime  is 
considered,  so  the  thrust  coeflBcient  and  forward  velocity  were  set  to  zero,  giving  an 
inflow  ratio  of  zero,  which  corresponds  to  a  steady  state  rotor  blade  pitch  angle  of 
(ft  )a  ~  0,  assuming  an  uncambered  blade.  This  simplifies  things  greatly  for  this  study  by 
ehminafrng  the  requirement  to  trim  the  rotor  S5^em  for  a  certain  aircraft  weight  and  flight 
condition.  The  aerodynamics  came  into  play  in  analyzing  the  effects  active  pitch  inpiits 
about  a  flat  pitch  condition  on  groimd  resonance  stability.  This  will  be  discussed  further  in 
a  later  section. 
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vm.  MOVING  BLOCK  TECHNIQUE 
One  of  the  drawbadcs  of  performing  direct  numerical  simulation  of  dynamic 
systems  is  that  time  histories  of  system  d^rees  of  freedom  only  offer  qualitative 
information  on  the  effect  that  certain  system  parameters  have  on  system  stability  or 
performance.  In  order  to  quantify  the  effects  of  varying  certain  system  parameters,  such 
as  rotor  speed,  flex-beam  stiffiiess,  and  active  control  inputs,  on  rotor-fiiselage  stability  in 
the  ground  resonance  regime,  a  method  was  needed  to  estimate  ^st^  damping  levels 
from  the  system  time  histories.  Moving  Block  Analysis,  a  technique  developed  at 
Locldieed  in  the  1970’s,  is  a  digital  method  of  analyzing  a  transient  time  history  to  obtain 
modal  damping  and  frequ^cy.  The  technique  is  first  described  in  some  detail  by 
Hammond  and  Dogget  [Ref  18]. 

The  technique  is  analytically  based  on  the  typical  transient  response  of  a  second 
order  systenL  Consider  the  following  transient  time  history, 

y(t)  =  A  sm(fi>r  +  &)  (8.1) 

where 

(o^  =G>„^(l-g^)  (8.2) 

The  finite  Fourier  trarrsform  of  XO  ^  the  danced  frequaaq?,  a ,  fix>m  time  r  to  r+  Tis, 

F(a>)  =  sin(mr  +  «9)  (8.3) 

After  carrying  out  the  integration  aiKl  making  the  assumptions  that  1,  so  ,  the 

magnitude  of  the  Fourier  transform  can  be  written  in  the  following  form. 
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where 


/(f)  =  +  (l  -  sm(2(®r + 5  ))  ■ 

e-'*"(l-c-‘'")fsiii[2(<o{r+r)+5)] 

Taking  the  natural  logarithm  of  equation  (8.4)  yields 


In|F(fi>)|  =  -gm + +  -Inf 
\l(oJ  2  \ 


The  last  term  in  equation  (8.6)  can  be  expanded  in  a  Taylor  series  to  yield 


ln|F(G>)|  =  -gm  + 

+ +  (e>7){sin(2(©r  +  5  ))  -  sin(2[0(T  +  7)  +  5  ])}] 

_  £  ^ j  sin(2(n)r  +  .9  ))  -  3  sin[2(Q>(r  +  T)  +  5  )] 

4  a)T+ sin(2{coT  +  5  ))  -  sin[2(6j(T  +  7)  +  5  )] 


From  equation  (8.7)  it  can  be  seen  that  if  lnjF(G))|  is  plotted  versus  r,  the  resulting  curve 
win  be  the  superposition  of  a  straight  line  with  slope  -  g-©  and  an  oscillatory  component 

which  osciUates  about  the  straight  line  with  frequency  2<y  .  If  h  is  assumed  that  Tis  an 
integral  multiple  of  the  basic  period  of  osciUation,  such  that 


27rN 


equation  (8.7)  reduces  to 


lnlF(ni)l  =  -gojT  +^$-sin(2(fiir  +  ,9)  +  c) 
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where  C  is  a  constant  given  by 


C  =  (8.10) 

From  equations  (8.9)  and  (8. 10)  it  can  be  seen  that  if  successive  discrete  Fourier 
transforms  at  a  frequency  &  are  paformed  for  0<  r  <  -  J ,  v^drere  r,  is  the  total  signal 

length,  a  plot  can  be  made  from  '\^ch  the  dan:q)ing  can  be  determined.  It  is  this 
procedure  ^ch  is  the  baas  for  the  moving  block  analyas  [Ref.  18]. 

For  a  sampled  signal,  the  moving  block  method  is  ^plied  by  first  estimating  the 
frequency  of  interest  embedded  in  the  signal  using  a  Fast  Fourier  Transform  (EFT).  A 
blodc  length  is  selected  consisting  of  Nb  data  points,  and  the  moving  block  fimction, 

ln|F(fi>,r^  ,  is  calculated  for  r  =  0 .  The  block  is  then  shifted  one  data  point  (time  step)  at 
a  time  and  the  moving  block  fimction  reconqjuted  for  r  =  » A/ ,  where  n  =  0, 1, 2, . . . ,  AT- 

Nb.  The  plot  of  lnjF(G>,  r)|  versus  t  is  fitted  with  a  linear  least  squares  fit,  and  tiie 
dancing  is  estimated  from  the  slope  of  the  curve  [Ref.  19]. 

For  this  study,  a  moving  block  analyas  program  was  developed  with  MATLAB®. 
The  m-file  code  is  contained  in  .^)pendix  E.  Moving  block  was  applied  in  the  code  by 
adapting  and  combining  the  procedures  outlined  by  Hammond  and  Doggett  [Ref  18]  and 
Bousman  and  \^^nkler  [Ref.  19].  The  method  used  for  computing  the  moving  block 
function  once  the  frequency  and  block  size  were  determined  was  to  evaluate  the  Fourier 
coefficients  for  the  first  time  block  with  the  following  relations 

«*(«■)  =  -^Z/(^+^)cos——  (8.11) 

^'b  x=0 


53 


(8.12) 


where 


6t(r)  = 


N,-\ 

x=0 


sin 


iTrk^x 

N, 


=a)Ni,M 


(8.13) 


zndfix^f)  is  the  signal  data  (for  the  first  block  r=  0) .  The  Fourier  coefficients  at  the 
next  time  step  are  then  calculated  using  the  followng  recur^on  relations, 

r 

2 


a^{z  + 1)  =  |at(r)+^[/(A^j  +  f)-/(r)]|cos^^^^' 


N. 


6,(r +  1)  =  -|a,(r)  +  +  ,)-/(^)]jsin[^]  +*,(r) 


1  I.  /  \  •  r 

J+».Ws.n[  J 

(8.14) 

J  \  Ny  J 

(8.15) 

The  magnitude  of  the  natural  logarithm  of  the  moving  block  fimction  is  then  given  by 


ln[F(fi>,r)]  =  jln[a*(r)'  +^(r)']  (8.16) 

The  accuracy  and  speed  of  the  moving  block  analysis  code  is  dependent  on  several 
fectors.  The  fi^uency  resolution  of  the  EFT  algorithm  is  inversely  proportional  to  the 
agnal  length  and  sampling  rate,  or. 


Am  = 


1 

NM 


(8.17) 


For  the  current  study,  rignal  fi^equencies  will  range  fix)m  approximately  5  to  30  Hz,  and 
sampling  rates  will  be  between  1 00  and  2000  Hz  with  the  typical  record  lengths  of  2  to  5 
seconds.  The  worst  resolution  considering  these  figures  would  be  approximately  0.5  Hz, 
an  error  of  1 0%  for  the  low  fi’equency  signal.  To  compensate  for  this,  a  refinement 
procedure  [Ref  18]  is  incorporated  into  the  code.  Additionally,  the  FFT  algorithm  is 
optimized  to  operate  on  record  lengths  in  powers  of  two.  If  the  record  length  is  not  a 
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power  of  two,  it  is  automatically  padded  with  the  required  number  of  zm>s,  w^di 
d^rades  the  accuracy  of  the  frequency  estimation.  To  remedy  this,  all  dgnal  lengths  were 
controlled  to  be  exactly  in  powers  of  two  by  adjusting  stop  time  and  the  ^e  of  time  steps 
when  ^ecuting  simulations. 

The  moving  block  code  developed  for  this  study  was  specialized  for  handling  imi- 
modal  or  bi-modal  agnals,  i.e.,  signals  with  one  to  two  dominant  modes.  For  the  bi-modal 
case,  suffidCTt  frequency  separ^on  must  exist  such  that  the  resolution  offered  by  the 
methods  described  above  will  be  adequate  ^ough  for  accurate  damping  estimates  of  both 
modes.  A  test  signal  of  the  following  form  was  used  to  verify  the  accuracy  of  program. 


/{>)  =  4e!cp 


(8.18) 


Parameter  values  were  fixed  to  match  the  test  case  considered  in  Ref  9,  and  were:  Ai = A2 
=  1000,  =  0.01,  $”2  =  0.02,  ©j  =  8.0  Hz,  and  ©2  =  6.0  Hz.  The  results  of  the  moving 

block  analyds  on  the  test  dgnal  shown  in  Figure  8. 1  is  shown  in  Iigure  8.2.  The  first  plot 
of  Figure  8.2  diows  the  resulting  power  spectrum  generated  by  a  FFT  of  the  test  signal. 
The  second  plot  is  die  result  of  refining  the  FFT  frequency  estimates.  The  third  plot  is  the 
crux  of  the  moving  block  analysis  vdiere  the  negative  of  the  slopes  of  the  straight  line  least 
square  fits  give  the  damping  ratios  f^en  divided  by  the  corresponding  damped 
frequencies  of  the  particular  mode.  The  damping  and  frequences  obtained  from  Figure 
8.2  are  summarized  in  the  Table  8.1. 
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Figure  8.1  Moving  Block  Test  Signal 


Figure  8.2  Result  of  Moving  Block  on  Test  Signal 
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Table  8.1  SamiDary  Results  for  Moving  Block  Bi-modal  Signal  Test  Case 


Parameter  -> 

^2 

Testdgnal 

8Hz 

0.01 

6Hz 

0.02 

Moving  Block  Analyds 

8.0314  Hz 

0.0098 

5.9S35  Hz 

0.0198 

From  the  results,  it  can  be  seen  that  the  two  modes  of  the  test  dgnal  are 
suffidendy  &r  apart  in  order  to  obtain  reasonably  accurate  dancing  estimates.  It  is, 
however,  in^mrtant  to  note  that  as  the  frequency  separation  of  a  bi-modal  dgnal  decreases 
to  approximately  5%  [Ref  19],  accurate  Hamping  estimates  win  no  longer  be  posdble. 

Figure  8.3  shows  the  effect  of  varying  rotor  speed  on  first  lead-lag  mode  damping 
as  determined  by  a  moving  block  analyds.  The  configuration  used  for  Figure  8.3  is  the 
basic  configuration  with  a  moderate  amount  of  damping  added  to  the  fiiselage  and  blade 
motions  in  order  keq)  time  histories  widiin  reasonable  bounds  when  simulating  indde  the 
self-exdted  r^on.  The  initial  exdtation  was  provided  by  setting  an  initial  fiiselage 
displacement  in  the  x-direction.  The  rotor  fiequenc^  (absdssa  in  Figure  8.3)  is  non- 


dimendonalized  by  the  fiiselage  natural  frequency  (,  ).  The  centCT  of  the  regresdng 


lead-lag  mode  instability  should  occur  at  a  non-dhnendonal  frequency 

_a  1 

Ho 


(8.19) 


for  a  rotor  ^stem  modeled  with  point  masses.  For  the  bade  configuration,  the  crater  of 

a 

instability  corresponds  to  -p—  =  1288 ,  which  is  in  good  agreement  Figure  8.3. 
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First  Lead-Lag  Mode  Damping  vs.  RotorSpeed 


Figure  8.3  Effect  of  Varying  Rotor  Speed  on  First  Lead-iag  Mode  Damping 

Figure  8.4  shows  the  results  of  a  moving  block  analysis  completed  for  rotor 
systems  with  different  Deutsch  numbers.  Deutsch  [Ref  21]  derives  a  criteria  to  determine 
the  quantity  of  damping  necessary  to  eKminate  ground  resonance  instability  through  the 
ftdl  range  of  rotor  speeds.  The  criteria  requires  that  the  product  of  the  blade  and  fuselage 
damping  parameters  be  greater  than  a  certain  parameter  determined  by  the  rotor-fuselage 
configuration.  The  rotor-fuselage  configuration  parameters  for  the  case  of  the  simple 
model  with  an  isotropic  pylon  and  rotor  are  given  by  the  following. 


(8.20) 


(8.21) 


(8.22) 
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(8.23) 


I  =  m,R^ 


The  dflTnping  parameters  for  the  p5ion  and  rotor  blades  are 

C„ 


^  {Mp  +  NiHj,)^ 


(8.24) 


1 - - 


(8.25) 


t^ere 


CD  = 


Doitsch’s  criteria  for  elimination  of  ground  resonance  is 


X  pX  A 

D=  .  "->1 


where 


P  = 


l  +  .^Aj  +  A2 

tta; 


(8.26) 


(8.27) 


(8.28) 


Here,  the  ratio  of  the  dancing  product  to  the  configuration  parameter  is  defined  as  the 
Deutsch  number ,  D . 
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Moving  Block  Results  Parametaed  With  Deulsch  Criteria 


Figure  8.4  Moving  Block  Results  Parametized  by  Deutsch  Number 
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IX.  ACTIVE  CONTROL 


Several  studies  have  addressed  the  problem  of  eliminating  helicopter 
af^otnpirbanir^l  instabilities  ^vith  active  control.  Straub  [Ref.  12]  and  Straub  and 
Wannbrodt  [Ref  13]  use  linear  state  space  methods  to  investigate  the  effects  of 
systematically  varying  feedback  gain  and  phase  in  a  closed  loop  ^stem  on  the  rotor- 
fuselage  dynamic  behavior.  Takahashi  and  Friedmann  [Ref  22]  move  one  stqp  fiirther  by 
proposing  a  simple  closed  loop  controller  based  on  an  optimal  state  estimator  in 
cor^ction  with  optimal  state  feedback  determined  from  linear  quadratic  Gaussaan  (LQG) 
optimization  techniques.  Weller  [Ref.  23]  ^owed  by  experiment  that  a  fixed  gain 
controller  which  transforms  fiiselage  states  into  swashplate  inputs  can  greatly  iiiq)rove 
asromechanical  stability  margins  and  eliminate  unstable  envelopes.  Wood ,  et  al.,  [Ref. 

24]  detailed  the  deagn  and  implementation  of  a  hi^er  harmonic  pitch  control  system  and 
demonstrated  that  it  can  be  an  effective  means  of  vibration  reduction.  For  the  scope  of 
this  study  a  similar  ^)proach  to  the  one  used  by  Wdler  [Ref.  23]  in  his  experimental 
investigation  was  incorporated  into  the  Emulation  enviromnent. 

In  the  most  general  case,  the  complex  modd  allows  rotor  blade  pitch  iiq)uts  to  be 
of  one  another  so  that  tiie  amulation  of  individual  blade  control  is  possible. 
The  general  fr>rm  of  a  pitch  input  b, 

=(^,)tSin(nfir+«>*)+(^^)^  cos(/int+€>t)  (9.1) 

where  w  is  the  harmonic  number  of  the  pitch  frequency,  nnd  the  iiq>ut 

phase  and  amplitude  weighting  parameters  (for  a  first  harmonic  input  they  would 
correspond  to  lon^tudinal  and  lateral  cyclic  fisr  the  k***  rotor  blade),  and  is  the 
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azimuthal  phase  angle  of  the  k*  rotor  blade.  For  the  following  active  control  example, 
active  pitch  inputs  via  a  swashplate  were  simulated  for  the  a  three  bladed  rotor,  thus,  n 
was  set  equal  to  unity  and  the  amplitude  and  phase  weighting  parameters  for  each  blade 
were  set  equal  to  each  other  such  that. 


(9.2) 

(9.3) 

Pylon  pitch  and  roll  position  feedback  was  transformed  into  swashplate  control 
iiq>uts  by  the  following  fixed  gain  relationship 


'g: 

'  Kcoi(f>)  ^sin(^)Tr,' 

A. 

-K^[<f>)  Kcoi^)\rt. 

(9.4) 


Stability  measurements  were  made  based  on  time  histories  of  the  orthogonal  components 
of  the  rotor  center  of  gra\aty  offset  given  by 


SK). 

^  =  1 


(9.5) 


y  eg 


*=1  * 


Zk) 

ir=1 


k 


(9.6) 


wliere, 

K)*  =  +co^V'*)«  (9'7) 

+  sin(^,)co^^JX^‘s)i  +  sin(v^*)e  (9.8) 
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These  time  histories  contain  both  r^essing  and  progresdng  mode  contributions.  The 
Hamping  levels  of  these  modes  for  various  gain  and  phase  settings  were  determined  using 
the  moving  block  analysis  program  described  in  a  previous  section. 

Initially,  a  simulation  was  run  with  the  feedback  gain  and  phase  set  to  zero  in  order 
to  get  a  baseline  response.  Figure  9.1  shows  the  results  of  the  basdine  response  of  the 
rotor  center  of  gravity  offeet  position.  Figure  9.2  displays  the  results  of  performing  a 
moving  blodc  analysis  on  the  rignal  of  Figure  9.1.  It  can  be  seen  in  both  figures  that 
the  center  gravity  of^  signal  contains  both  the  r^resring  and  progressing  lead-lag 
modes.  The  progressing  mode  danq>s  out  relatively  quicldy  and  is  of  little  influence  after 
approximately  0.4  seconds  of  simulation  (h  can  be  seen  that  the  high  fi'equ^cy  con[q}onent 
of  the  rignals  in  Figure  9.1  “washes  out”  by  0.4  seconds). 
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rad(aRsA2/Hz 


figure  9.1  Rotor  c.g.  Offset  For  Baseline  Case  (K=0,  =0) 


figure  9 J  Moving  Block  Results  for  Baseline  Simulation  (K=0,  (j)=0) 
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The  moving  block  analysis  determined  a  damping  ratio  of -0.0247  (©-7.1511  Hz) 
for  the  regressing  mode  and  a  damping  ratio  of  0.0230  (©=27.1916  Hz)  for  the 
progressing  mode.  The  next  step  was  to  conduct  the  same  analysis  as  was  completed  on 
the  baseline  case  for  a  range  of  active  control  phase  angles  at  fixed  values  of  gain. 

Figure  9.3  displays  the  results  of  running  a  controller  phase  sweep  at  gain  settings 
of  K=  0.4  and  Ar=0.8.  The  case  chosen  is  slightly  unstable,  with  a  no  control  damping 
ratio  shown  by  the  green  line.  The  results  demonstrate  that  stability  can  be  improved  by 
active  swashplate  control  inputs  and  that  the  simulation  techniques  used  in  this  study  can 
provide  a  useful  tool  for  predicting  which  gain  and  phase  combinations  would  be  required 


for  a  single  controller. 


Itgare  9.3  Damping  Ratio  For  Controller  Phase  Sweep  (K“0.4  and  K=0.8) 

Figure  9.4  shows  the  SMULINK®  model  utilized  for  complex  model  simulations 


and  the  active  control  analysis. 
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Figure  9A  SIMULINK®  model  implementing  fixed  feedback  gain  active  rotor  control 
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X.  CONCLUSIONS  AITOWECOMMENDATIONS  FOR  FUTDEE  RESEARCH 
A.  CONCLUDING  REMARKS 

A  method  for  fonmilating  and  automatically  coding  the  equations  of  motion  of  a 
coupled  rotor-fiiselage  system  1^  use  of  symbolic  processing  software  and  (fynamic 
gtniilatinn  software  has  been  developed.  The  resulting  mathematical  models  were  used  to 
perform  amulations  of  coupled  rotor-fiiselage  ^sterns  in  ground  resonance.  Analysis  of 
the  Hynamic  and  stattility  diaracteristics  were  quantified  using  the  moving  block  tedinique. 
A  simple  rotor  modd  was  used  to  dmnonstrate  essential  charactoistics  of  ground 
resonance  and  the  effects  that  parameter  variations  such  as  rotor  speed,  flekbeam  dastic 
behavior,  damp^  feilure,  and  rotor  blade  damage  have  on  those  characteristics.  A  more 
complex  model,  adding  fuselage  phdi  and  roll  and  rotor  blade  flap  degrees  of  fi-eedom, 
was  used  to  demonstrate  how  the  modeling  tedmique  could  be  used  to  mqilore  the  effect 
of  active  rotor  control  on  ground  resonance.  The  modeling  technique  proved  to  be  a  very 
powerful  tool  in  that  it  eliminated  the  time  consuimng  process  of  manualfy  driving  and 
coding  the  very  complex  equations  of  motion  of  a  multi-degree  of  fi'eedom  rotor  ^stem 
into  a  dynamic  simulation  mivironment.  By  integrating  SIMULINK®  into  the  process, 
with  its  versatility  in  analysing  dynamic  systems,  the  technique  has  direct  application  to  the 
dedgn  of  modem  damperless  rotor  systems. 
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B.  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 


1 .  Addition  of  rotor  blade  torsional  degrees  of  freedom  to  overall  model. 

2.  Addition  of  geometric  characteristics  such  as  pre-cone,  pre-sweep,  ofi^  hinge 
inclination,  elastic  coupling  (pitch-lag,  lag-flap,  etc.). 

3 .  Addition  of  unsteady  aerodynamics  sudi  as  a  finite  state  wake  model  or 
dynamic  inflow  model. 

4.  Incorporation  of  a  trim  routine  so  the  model  can  be  used  for  hover  and  forward 
flight  aeromechanical  stability  analysis. 

5.  A  comprehensive  study  of  active  control  methods  using  the  developed 
modeling  technique,  including  optimization  techniques  (such  as  LQR  and 
LQG)  and  individual  blade  control 

6.  Validation  of  simulation  results  with  experiment. 
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APPENDIX  A 


MAPLE®  WORKSHEET  USED  TO  GENERATE  AND  CODE  TBffi  EQUATIONS 
OF  MOTION  FOR  A  COUPLED  ROTOR  FUSELAGE  SYSTEM 
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APPENDIX  A 


(MAPLE  WORKSHEET) 

EQUATIONS  OF  MOTION  FOR  A  HELICOPTER  IN  GROUND  RESONANCE 
CONSIDERING  4  FUSELAGE  DEGREES  OF  FREEDOM  AND  ROTOR  BLADE  FLAP  AND 

LEAD-LAG  DEGREES  OF  FREEDOM 

[  >  restart: 

'  >  withdinalg)  : 

Warning,  new  definition  for  norm 
.  Warning,  new  definition  for  trace 
[  >  diffl:  =  (arg)->map(diff ,arg,t) : 

[  Define  coordinate  transformations: 

[ - - - 

>  psi :=Omega*t+Phi [k] ; 

\\f  + 

>  T3:*=alpha->matrix(3,3,  [cos  (alpha)  , sin  (alpha)  , 0 , -sin  (alpha)  , cos  (alp 
ha) ,0,0,0,!])  ; 

T3  :=a-^  matrix(3,  3,  [cos(aX  sin(a),  0,  -stn(a),  cos(a),  0,  0,  0, 1  ]) 

>  T2  :=alpha->inatrix (3,3,  [cos  (alpha)  , 0,  sin  (alpha)  ,0,1,0,  -sin  (alpha)  , 0 
, cos (alpha) ]) ; 

72  :=  a  matrix(3,  3,  [cos(a),  0,  sin(a),  0,  1,  0,  -sin(aX  0,  cos(a)]) 

>  Tl:=alpha->matrix(3,3, [1,0, 0,0, cos  (alpha) , sin (alpha) , 0 , -sin (alpha) 
,cos (alpha) ] ) ; 

7i  :=  a  ^  matrix(3,  3,  [  1,  0,  0, 0,  cos(a),  sin(a),  0,  -sin(a),  cos(a)]) 

>  Ml:=transpose(multiply(Tl(r[l] (t) ) ,T2(r[2] (t) ) ) ) ; 

'cos(r^{t))  -sin(r,(0)sin(r2(0)  -cos(ri(0)  sin(r2(0) 

MI  :=  0  cos(7'i(0)  -sin(ri(0) 

_sin(r2(0)  sin(r,(0)cos(r2(0)  cos(ri(/))  cos(  7*2(0 ) 

>  M2 :=transpose (T3 (psi)) ; 

cos(%l)  -sm(%l)  0 

M2  :=  sin(%l)  cos(%l)  0 

0  0  1. 

%i  :=  a  /  + 

>  M3 :=transpose (multiply (T3 (zeta [k] (t) ) , T2 (beta [k] (t) ) ) )  ; 

cos(C,(0)cos(|3^(0)  -sin(Cfc(0)cos(|3,(0)  -sin(p/0) 
sin(Ct(/))  cos(Cfc(0)  0 

cos(Cfc(0)sin(pfc(0)  -sin(Cfc(0)sin(P;k(0)  cos(Pfc(0) 
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M3:= 


■  >  M4:=multiply(T3(zeta[k]  (t) ) ,T2 (beta[k]  (t)) ,T3(psi) ,Tl(r[l]  (t) ) ,T2( 
r[2] (t))) : 

[ - 

[  Energy  expressions  for  kth  rotor  blade 

[ - 

[  Kinetic  energy  of  kth  blade  (TBk) 

[ _ 

[  >  rhoFI__I :  =vector  ( [u[l]  (t)  ,u[2]  (t)  ,0] )  : 

[  >  rhoHF : =vector ( [ 0 , 0 , h ] ) : 

[>  rhoHF_I:==multiplY{Ml,rhoHF)  ; 

[  >  rhoBuH : ^vector ( [ el , 0 , 0 ] ) : 

[  >  rhoBuH_I :  =multiply  (Ml  ,M2  ,  rhoBuH)  : 

[  >  rhoPBd : =veator ( [R, 0 , 0 ] ) ; 

[  >  rhoPBd_I : =mul tiply (Ml , M2 , M3 , rhoPBd) : 

[  >  rho  :==matadd  (rhoFI_I  ,niatadd  (rhoHF_I  ,matadd  (rhoBuH_I ,  rhoPBd_I) ) )  : 

[  >  V:=diffl (rho) : 

[  >  Vsqr:=V[l]'"2+V[2]^2+V[3]''2: 

■>  TBk:=l/2*iiib[k]*Vsqr/ 

1  (((d  ^  (d  '1 

TBk  =-mbA  -«i(0  +sin(ri(0)  -^,(0  sinC/^CO)/* 

2  \\\ot  J  \Ot  j 

( d  \  (  ^5  ^ 

-  cos(ri(0)  cos(r2(0)  A  +  -sin(r2(0)  T^2(0  cos(%3)  -  cos(r2(0)  sin(%3)  n 

\(jt  J  \  \Ol  J 

(d  \  r 5  'I 

-  cos(r.(0)  -r^t)  sin(r2(0)  sin(%3)  -  sin(ri(0)  oosir^{t))  sin(%3) 

\dt  J  ^ 

-  sin(ri(0)  sin(r2(0)  cos(%3)  eJ  +  |^|^-sin(r2(/))  J  cos(%3) 

-  cos(r2(/))  sin(%3)  Q  -  cos(ri(0)  sin(r2(0)  sin(%3) 

fd  ^  \ 

-  sin(ri(0)  cos(r2(/))  -r2(0  sin(%3)  -  sin(ri(0)  sin(r2(0)  cos(%3)  Q J  cos(C*(0) 

\dt  J 

cos(P;i.(0)  -  (cos(r2(/))  cos(%3)  -  sin(A'j(0)  sin(r2(/))  sin(%3))  siii(C^(0)  %2  cos(3*(0) 

f 

-  (cos(r2(0)  cos(%3)  -  sm(ri(0)  sm(r2(0)  sin(%3))  cosCCfcCO)  sin(p;.(/))  %1  +  ^ 

( d 

sin(/2(0)  T''2(0  sin(%3)  -  005(7^2(0)  cos(%3)  Q 
\dt  J 

(d  \  1 

-  cos(7‘i(0)  J  sin( 7-2(0)  cos(%3)  -  sin(7-i(0)  cos(r2(0)  J  cos(%3) 

■n 

+  sin(7-j(0)  810(7-2(0)  sin(%3)  fl  sin(C*(0) 
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+  (-cos(r2(/))  sin(%3 )  -  sin(r,(0)  stn(r2(0)  cos(%3))  cos(Cfc(0) 

f  d 

+  sin(ri(0)  sm(r^(t))  cos(Q.(0)  sin(P/XO) 

\dt  J 

( d  ^ 

-  cos(ri(0)  cos(r2(/))  cos((;,XO)  sin(P^(0) 

+  cos(ri(/))  sin(r2(0)  sin(Q(0)  sin(p^(0) 

1  Y  ffd  1  fd  ) 

-cos(/-i(0)sin(r2(/))cos(CA:(0)cos(|ijt(0)%lJ-^j  j-cos(ri(0)[^^^^i(0 

-sin(rj(/))  sin(%3)e7  +  cos(rj(0) cos(%3) D + 

fd  ^ 

-sm(r,(0)  T^i(0  sin(%3)  cos(C;t(0)  cos(P;t(/)) 

\dt  J 

+  cos(ri(0)  cos(%3)  Q  cos(Ct(0)  cos(p^(0)  -  cos(r,(0)  sin(%3 )  sin(C/0)  %2  cos(Pfc(0) 

(d  ^ 

-  cosirXO)  sm(%3)  cos(Cfc(0)  sin((5fc(0)  %1  “  sin(rj(0)  J  cos(%3)  sin(C^(0)  , 

-  cos(rj(0)  sin(%3)  n  sin(C*(0)  +  cos(ri(0)  cos(%3)  cos(Cjfc(0)  “/»2 

+  sin(Cjt(0)  %2  sin(p/XO)  sin(r,(0)  -  cos(Cfc(0)  cosO^CO)  %1  sin(r,(/)) 

fd  y  f  fd  ) 

-cos(Cfc(0)sin(P/0)cos(ri(0)[^“^i(0 +  |^-sin(ri(0)  jcos(72(0)^ 

f  d  ^  f  fd  \ 

-  cos(ri(0)  sin( 7-2(0)  [^O(0j  h  +  j^cosC/jCO)  [^O(0 J  cos(%3)  -  sin(r2(0)  sin(%3)  O 

f d  ^  ^  . 

+  cos(ri(/))  cos(/-2(/))  sm(%3)  -  sin(rj(/))  sin(r2(0)  J  sin(%3) 

^  rr  1 

+  sin(ri(0)  cos(r2(0)  cos(%3)  QJ  el  +  |^[^cos(r2(/))  (^^0(0 J  cos(%3 ) 

f  d  \ 

-  sin(r2(/))  sin(%3)  Q  +  cosir^t))  [^O(0j  oos(r^{t))  sm(%3) 

-  sin(ri(0)  sin(r2(/))  rO(0  sin(%3 )  +  sin(7-,(/))  cos(r2(0)  cos(%3 )  Q J  cosC^^-CO) 

cos(p^(/))  -  (sin(/-2(0)  cos(%3 )  +  sin(/-,(0)  cos(r2(/))  sin(%3))  sin(^^(0)  %2  cos(P^(0) 

/ 

-  (sin( 7-2(0)  cos(%3 )  +  sin(7-i(/))  cos(r2(0)  sin(%3))  cos(Ci(0)  sin(pfc(0)  °/«l  +  ^ 
-cos(7'2(/))  I  ^0(Oj  sin(%3)  -  sin( 7-2(0)  cos(%3)  O 
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fs  1  1 

+  cos(ri(0)  cos(r2(0)  cos(%3)  -  sin(rj(0)  sin(72(0)  o/2(^)  I  cos(%3) 

\ut  J  \0l  J 

\ 

-  sin(/'i(0)  cos(r2(0)  sin(%3)  sin(Cjt(0) 

+  (-sin(r2(0)  sin(%3)  +  sin(rj(0)  cos(7'2(0)  cos(%3))  cos(C;t(0)  %2 

-sin(/-i(/))  “/'i(0  cos(r2(0)cos(Cfc(0)sm(p^(0) 

\ot  j 


r 


d 


-  cos(ri(/))  sm(/:2(0)  J  cos(C*(0)  siti(P*(0) 

-  cos(ri(?))  cos(72(/))  sm(C*(0)  %2  sin( (3^.(0) 


^2^ 


+  COS(7-i(0)  COS(7'2(0)  COS(Ca:(0)  COS(|3^(0)  %1  U? 


%2:=|«0 
%3  :=  Q  ^ 

c _ 

[  Potential  energy  of  kth  blade  (UBk) 

[ _ 

[  >  UBkl:=l/2*(beta[k]  (t) ''2*kfl  [k] +zeta[k]  (t) ''2*kll  [k] )  : 

[>  UBk3 :  =1/4 *  (beta [k]  (t)  ^4*kf3  [k] +zeta(k]  (t)  ^4*kl3  [k] )  : 

■  >  UBk : =UBkl+UBk3 ; 

UBk  +  \  WO*  + j«0*  kih 

c _ 

[  Dissapative  energy  of  kth  blade  (pBk) 

[ _ 

r  >  DBk :  =1/2 *diff  (beta [k]  (t)  , t) '''2*cf  [k] +l/2*diff  (zeta [k]  (t)  ,t)'^2*cl[k] 


DBk=- 


1 


^2 


\  d 


cL 


[ - 

[  Energy  expressions  for  fuselage 

[ - - - 

[  Kinetic  energy  of  fuselage  (TF) 

[ _ 

[>  TFt:=l/2*  (diff  (u[l]  (t)  ,t)'"2*M[l]+diff  (u[2]  (t)  ,t)^2*M[2])  : 

[  >  TFr:=l/2*  (diff  (r  [1]  (t)  ,  t)  ^2*I11+H^^f^^^(r  [2]  (t)  ,t)  ^2*I22-2*diff  (r[ll  ( 


[  t)  ,t)*diff  (r[2]  (t)  ,t)*I12)  : 
>  TF:=TFt+TFr; 


1 

7F:=- 

2 


+o 


f  . 


^2 


f  . 


112 


1 

M,+- 


1U.\ 


^2 


^07 


722 


[ _ 

[  Potential  energy  of  fuselage  (UF) 

[ _ 

[  >  UFt:=l/2*u[l]  (t)''2*KTl+l/2*u[2]  (t)'^2*KT2; 
[  >  UFr:=l/2*r[l]  (t) ^2*KRl+l/2*r [2 ] (t)^2*KR2: 
r  >  UF:=UFt+UFr; 


VF :=^u,( t)  KTl  +  ^ tf  KT2  +  ^ KRl +\ r^{tf  KR2 

[ _ 

[  Dissapation  energy  of  fUselage  (DF) 

[  >  DFtv:=l/2*diff  (u[l]  (t)  ,  t) ‘^2*CTl+l/2*diff  (u  [2  ]  (t)  ,t)^2*CT2: 

[  >  DFrv:=l/2*diff  (r  [1]  (t)  ,  t)  ^2*CRl+l/2*diff  (r  [2]  (t)  ,t)''2*CR2: 

>  DFth:=l/2*diff  (u[l]  (t)  ,t)'^2*abs(diff  (u[l]  (t)  ,  t) )  *VTl+l/2*diff  (u  [2] 
(t)  ,t)'"2*abs(diff  (u[2]  (t)  ,t)  )*VT2: 

>  DFrh:=l/2*diff  (r  [1]  (t)  ,  t)  ''2*abs  (diff  (r[l]  (t)  ,  t)  )  *VRl+l/2*diff  (r  [2] 
(t) ,t)^2*abs(diff (r[2] (t) , t) ) *VR2 : 

>  DF : =DFtv+DFrv+DFth+DFrh ; 


l(d 

DF:=-\-um 

iKdt  *  J 

TWi(0 
\dt 


^2  \(  ri  V  >  1 

cr7+- 


_a 

2\.dt 


«2(0 


CT2+- 


2\.dt 


h{t) 


CR1+- 


1 

2 


d 


dt 


1  (  d 
VT1+- 

2 


VR2 


d 

dt 


VT2+- 


1 


(d  '] 

2 

CR2 

(d  > 

d  ,  , 

VRl 


[  Aerodynamics  (Generalized  Aerodynamic  Forces) 

[ 


>  Vair  :=vector  ( [niu*Omega*R,  0  ,Omega*lambda*R] )  ; 


L  Ffif/r  ;=  [p  Q7?,  0,  Q  Xi?] 

[>  V_I_t:=matadd(-V,Vair)  : 

[>  V_Bd_t :=map (simplify, mul tiply (M4, V_I_t) )  : 
[ >  UR:=V_Bd_t[l]  : 

[  >  UT:=-V_Bd_t[2] : 

[  >  UP:=-V_Bd_t[3] : 

[>  UU:=sqrt(UP'"2+UT'^2)  : 
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[  >  aoa:=theta[k] -arctan(UP/UT) : 

[  >  dFbeta:=l/2*rhol*a*c* (aoa*UU*UT-cdO/a*UU*UP) : 
[>  dFzeta:=-l/2*rhol*a*c* (aoa*UU*UP+cdO/a*UU*UT) : 


>  Mbeta_k:=0.7*R''2*dFbeta*cos  (zeta[k]  (t))  / 
Mbeta  k  :=  .3500000000 


K"  p\  a  c  -  arctan 
%1  Cl  t 


l%2^  +  %3^  %3 


cd0'J%2'^  +  %3^  %2 


cos(C*(0) 


%2  ~  R  cos(  C(  / ) )  T  W  0  +  cos(  r^(  t) )  sin(  /) )  sin(  %  1 )  sin(  r^{t))aXR 
\ot  J 


+  sin( p/ /) )  cos( %  1 )  sin( r^{t))Cl'KR-  sin( P*( 0 )  sin( %  1 )  cos( r-^{t))  ti^i t) 


-  cos(P,(/))  sin(ri(0)  -%(0  +  sinCr^CO)  sin(p,(/))  sin(%l )  sin(ri(0) 


-  sin(pi.(0)  cos(%l)  cos(/-2(0)  -«i(0  -  siti(/2(0)  cos(P^.(0)  cos(r,(0)  “  Wi(0 


-cos(r2(0)  cos(p^(0)  cos(ri(0)n  ?l;?  +  cos(3^(0)  co&{r^{t))el  I  cos(%l) 


+  cos( P^.( r) )  - ri( 0  sin(%l )  el  +  cos(  0 )  ^  ^ ^i( 0  cos(%l )  sin(Ci( /) ) 


+  sin(3^(0)  sin(%l )  -ri(0  U  +  sin(p/X0)  cos(%l )  cos(ri(0)  -^jCO  J  h 


+  sin( p^( /))  sin(ri(  t) )  R  sin(Ci.(  0 )  T  '’jC  0  +  cos(%l )  R  cos{r^{  0 )  “ r^i  0  j  cos(Cfc( 0) 


+  sin(%l  )i?  1^— ri(0j  cos(Cfc(0)  +  sin(3^(0)-^  sin(Q.(0)  ^ 

( d  '\ 

+  sin(p^(0)  cos(%l)  cos(r2(0)  \iClR-  cos(p^(0)  cos(ri(0)i?  sin(C^(0)  [^'^2(0 J  sin(%l) 
-  sin(r2(0)  sin(p/0)  sin(%l)  sin(ri(/))  |i  Hi?  +  cos(pfc(0)  cos(ri(0) 


%3  :=  sin(4(0)  cos(3^(0)  sin(%l)  ^-rjCOj  h  +  cos(C^(0)  elCl  +  co^{^J,t))RCl 


+  — Cfc(0  +sin(C*(0)cos(p^(/))cos(%l)cos(ri(0)  /2(0j^ 


sm(p/0)^  -/xU)  cos(%l)-sin(C^.(0)sin(p,(0)  -/'i(0  sin(%l)e7 


(d 

(d  \ 

-7*2(7) 
\dt  ’) 

+  cos(Q.(7))  sin(%l )  oos( 7*1(7)) 

\dt  J 

d  ]  f  d  ] 

-  cos(C^;(0)  cos(%l )  j  h  -  sin(C/XO)  sin( (3^(0)  cos(ri(/))  el  |^^^2(0 J  cos(%l ) 

fa  ^  (d  ' 

+  sin( p^( 0 )  cos( t\{t))R  sin( %  1 )  +  sin( /-j ( 0 )  cos( P^( / ) ) i?  r2( 0 


\d( 


Kdi 


f  d  ^ 

-  cos(r2(0)  cos(Cfc(0)  sin(%l)  |^-//i(/)J  +  cos(r2(0)  sin(4(0)  cos(p^(0)  cos(%l)  i? 

+  cos(r2(/))  cos((^(t))  sin(%l)  ^  07? 

-  sin(72(0)  s^T-iCO)  sin(C/0)  cos(p;t(0)  sin(%l)  07? 

+  sin(r2(/))  sin(/'j(0)  cos(C*(0)  cos(%l)  07? 

-  sin(r2(0)  sin(C/XO)  sin(Pi.(7))  cos{r^(t))  |.i  07? 

(d  '] 

-  008(7*2(0)  sin(C//))  cos(p^(7))  cos(%l )  [^Wj(/) J 


+  sin( 7*2(0)  sin(7*i(0)  sin(Ci.(0)  cos(p^^(0)  sin(%l )  u^(t) 


-  sin(7‘2(/))  8171(7*1(0)  COS(Cfc(/))  C0S(%1  )  1^—  77i(0^ 

(d 

+  sin( 7*2(7))  sin(C;t(0)  sin(Pi.(7))  cos(7*i(7)) 


^dt 


TJhiO 


-  cos( 7*1(7))  sm(Cfc(7))  cos(P,(7))  sin(%l) 

C  a  'A 

+  cos( 7*1(7))  cos(^.(7))  cos(%l )  ~  7<2(7)^  +  sin((^^(7))  sin(P//7))  sin( 7*1(7)) 


f- 

^dt 


TJhO) 


+  sin(r2(7))  sin((!^^(7))  cos(P/7))  cos(%l )  O  XT?  +  sin( 7*2(7))  cos(Ct(0)  sin(%l)  O  XT? 
+  008(7*2(7))  8in(ri(7))  sin(C^(7))  cos(P;i.(7))  sm(%l) O  XT? 

-  cos(r2(7))  sin( 7*1(7))  cos(Cfe(7))  cos(%l)  O  XT? 

+  008(7*2(7))  8in(C/X7))  sin(pjt(7))  cos(7*i(7))  O  XT? 

>  Mzeta_k :=0 . 7*R^2*dFzeta; 

Mzetaji  \— 

r 


-.3500000000  7?^  pi  «c 
%1  :=0  7  +  <J)* 


VV 


0^  -  arctani 


%3 


^%2j 


V%3^+%2^  %3  + 


c^qV%3^  +  %2^  %2 


o 


%2  :=  sin(C;i(0)  cos(P^(7))  sin(%l)  ^  +  oos(Ck(0)  O  +  cos(p^.(7))T?  O 
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fd  )  (^1 

+  i?  1^—  Ci.( 0 J  +  sin(C*( 0 )  cos( |3^.( 0 )  cos(%l )  cos(ri(i ) )  ^2(0  j  ^ 

f  d  ^  f  ^  ^ 

-sin(p^(/))i?  -ri(0  cos(%l)-sin(C^(0)sin(p//))  -ri(0  sin(%l)^y 

\Ot  J  \0l  J 

(d  \  ^ 

+  cos(C^(0)  el  sin(ri(0)  [^'2(0j  +  cos(C,(0)  sin(%l )  cos(ri(/))  J  ^ 

(d  \  1 

-  cos(C^(0)  C0S(%1 )  ^  -  sin(C/XO)  sin(p^(0)  cos(ri(0)  U J  ^ 

fd  ^  1 

+  sin(  p^(  t) )  cos( r,( 0 )  r^(  0 J  sin(%  1 )  +  sin(rj( / ) )  cos(  Pa:( 0 )  ^  ^2i  0^ 

f  d  ^ 

-  cos(r^(t))  cos(C/XO)  sin(%l )  |^-  «i(0j  +  cos(a'2(0)  sin(C/c(0)  cos( p^(0)  cos(%i )  n  i? 

+  cos(  r2(  / ) )  cos(  Cfe(  0 )  sin(  %  1 )  |i  O  7? 

-  sin(r2(0)  sin(ri(0)  sin(C^(0)  cos(P/0)  sm(%l )  p,  i? 

+  5111(^2(0)  sin(ri(0)  cos(Cj.(/))  cos(%l )  p  O  7? 

-  sin(/-2(/))  sin(Cjt(0)  sm(3*(7))  cos(ri(0)  P 

-cos(r^(t))  sinCC/XO)  cos(p^(0)  cos(%l) 

) 

+  sm(r2(0)  sin(ri(0)  sin(C/X0)  cos(  p^.(0)  sin(%l )  Wi(/)^ 

fd  ) 

-  sin(r2(0)  sin(ri(0)  cos(i!^//))  cos(%l )  Hi(0^ 

fd  ) 

+  sin(r2(7))  sin(C^(0)  sin(pi(0)  cos(ri(/)) 

-  cos(ri(0)  sin(Ct(0)  cos(  pfc(7))  sm(%l )  |^— 

f  d  ^  ^ 

+  cos(ri(0)  oos(Q(t))  cos(%l )  +  sin(^k(0)  sin(P*(0)  sin(rj(7)) 

+  sin(r2(7))  sin((I^^(?))  cos(P^(7))  cos(%l )  fi  A,7?  +  sin(r2(7))  cos((!^^(0)  sin(%l )  Q  XT? 

+  cos(^2(7))  sin(rj(7))  sin(Cfc(7))  cos(P;t(0)  sin(%l)n  XT? 

-  cos(r2(t))  sin(ri(0)  cos(C^(7))  cos(%l )  Q  XT? 

+  cos(r2(0)  sin(?^^.(0)  sin(P^(/))  cos(rj(0)  ^ 

fd  ^ 

%3  :=  T?  cos(^/t))  P^( 0 J  +  cos(/-2( t) )  sin( p^(/) )  sin(%l )  sin(ri( /) )  X T? 


sin(p^(/))  cos(%l)  sin(r2(0)  sin(p^(0)  sin(%l )  cos(/'i(0) 


) 

J 


(d  ^ 

^  d  '' 

“«2(0 

\dt  J 

+  sin(/'2(0)  sin(3/X0)  sin(%l )  sin(rj(0) 

Kdt  ’) 

-  sin(P^.(/))  cos(%l )  cos(r2(0)  |  ) "  sin(r2(0)  cos(i5^.(/))  cos(r,(/))  Wi(0 


) 

J 


-  005(^2(0)  cos(p^XO)  ^os(/’i(/))  QXR-\-  cos(p^(/))  cos(r  ^(t))  el  ^^^2(0  j  cos(%l) 


J 


+ 


cos(  15^(0)  [f  ^i(Oj  sin(%l )  el  +  cos(p*(0)^  cos(%l )  sin(C/XO) 


(d 


\ 


(d  ^  f  5 

+  sin(p^(0)  sin(%l )  -r^{t)  h  +  sm((3^(0)  cos(%l )  cos{r^{t))  “^2(0 


\dt 


(d  ^  ^  \ 

+  sin(P^.(0)  sin(ri(0)-R  sin(C;fe(/))  [^'’aCOj  +  cos(%l  )i?  cos(ri(0)  [^'*2(0 J  cosiQ(t)) 


+  sm(%l ) -ri(/)  cos(C*(0)  +  sin(Pfc(/))  R  sin(Q(0)  ^ 
\dt  J 


fd 


sin(p^(0)  cos(%l)  cos(r2(0)  \i^R-  cos(%(t))  cos{r^{t))R  sin(C/XO)  J  sin(%l) 


[  -  sin(/-2(0)  sin(p/0)  sin(%l )  sin(r,(0)  li^R  +  cos(pi.(0)  cos(/’i(0) 

[ - - - - - - - - 

[  Derivation  of  equations  of  motion  by  Lagrangian  method 

[  This  section  defines  vectors  of  displacement  degrees  of  freedom,  their  rates  and  accelerations. 
[  >  DOFF:  =  [u[l]  (t)  ,u[2]  (t)  ,r  [1]  (t)  ,r  [2]  (t)  ]  : 

'  >  N :  =3  /  Choose  number  of  rotor  blades 

N:=3 

[  >  DOFB :  =  [ ] : ThetaB :  =  [ ] : 

'  >  for  i  from  1  to  N  do 
>  DOFB := [op (DOFB) ,beta[i] (t) ,Z6ta[i] (t) ] : 

_  >  od: 

■>  DOF  :  =  [op  (DOFF)  ,  op  (DOFB)  ]  ; 

DOF  :=  2/2(0,  ^i(0,  /^(O,  3i(0,  Ci(0,  PaCO,  ^^(0,  CaCOl 

r  >  dDOF:=diffl (DOF)  / 


dDOF  = 


[  >  ddDOF:=diffl (dDOF) ; 
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ddDOF  := 


d 


dt^  ^  If  a?  If  a/^ 

r  This  section  defines  transformations  between  time  dependent  and  independent  notation  in  terms  of 
substitution  sets. 


[>  setA;={ } : setB:={ } : setC:={ } : 

[>  setD:={} :setE:={} :setF:={} : 

[  >  DOFq:  =  [] :dDOFq:  =  [] :ddDOFq:  =  [] : 

’  >  for  i  from  1  to  vectdim(DOF)  do 

>  DOFq:=[op(DOFq) ,q[i] ] : 

>  dDOFq:=[op(dDOFq) ,dq[i] ] : 

>  ddDOFq : = [ op (ddDOFq) , ddq [ i ] 3 : 

>  setA:=setA  union  {ddDOF [i]=ddDOFq [i] } : 

>  setB:=setB  union  {dDOF[i]=dDOFq[i] } : 

>  setC:=setC  union  {DOF [i ] =DOFq[i] } : 

>  setD:=setD  union  { ddDOFq [ i ] =ddDOF [ i ]} : 

>  setE:=setE  union  {dDOFq[i]=dDOF[i] } : 

>  setF:=setF  union  {DOFq[i]=DOF[i] } t 

_  >  od: 

r  >  setl:=setA  union  setB  union  setC ;  Substitution  set  to  go  fi’om  dependent  to 


L  ind^endent 

>  set2:=setD  union  setE  union  setF:  Substitution  set  to  go  firom  independent  to 
dependent 

[  This  section  combines  all  of  contributions  to  the  terms  of  the  Lagrange  equation. 


[  >  T;=TF: 

[  >  U;=UF: 

[  >  D1:=DF: 

[  >  GF:  =  [0, 0,0,0] : 
r  >  for  i  from  1  to  N  do 


>  T:=T+subs (k=i,TBk) : 

>  U:=U+subs (k=i,UBk) : 

>  D1 :=Dl+subs (k=i,DBk) : 

>  GF : = [op (GF) , subs (k=i ,Mbeta_k) , subs (k=i ,Mzeta_k) ] : 

_  >  od: 

[  This  section  carries  out  the  differentiation  operation  of  the  Lagrange  equation  one  term  at  a  time 
[>  Temp:=subs  (setl ,  T)  : 

■>  for  i  from  1  to  vectdim(DOF)  do 

>  tempi :=diff (Temp,dDOFq[i] ) : 

>  temp2 : =subs ( S6t2 , tempi ) : 

>  temp3:=diff (temp2, t) : 

>  LI :=subs (setl , temp3) : 

>  L2 :=diff (Temp,DOFq[i] ) : 
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>  L3:=diff (subs (setl,U) ,DOFq[i]) : 

>  L4 ;=diff (subs (setl,Dl) ,dDOFq[i] ) : 

>  GFq:=subs (setl,GF[i] ) : 

>  EOM[i] :=Ll-L2+L3+L4-GFq: 

_  >  od: 

[  This  section  formats  the  equations  of  motion  into  the  form  A  d2x/dt2  =  f 

'  >  A:=matrix  (vectdim(DOF)  ,vectdiin(DOF)  )  ; 

A  :=  array(l ..  10,  1  ..  10,  [  ]) 

■  >  for  i  from  1  to  vectdim(DOF)  do 

>  for  j  from  1  to  vectdim(DOF)  do 

>  A[i, j] :=coeff (EOM[i] ,ddDOFq[j]) : 

>  od: 

_  >  od: 

[  >  sotZ :  =  { } : 

'  >  for  i  from  1  to  vectdim (ddDOFq)  do 

>  setZ:=setZ  union  {ddDOFq[i]-0} : 

_  >  od: 

'  >  f :=array (1 . .vectdim (DOF) ) ; 

/:=array(l  ..  10,  [  ]) 

■  >  for  i  from  1  to  vectdim (DOF)  do 

>  f [i] :=-eval (subs (setZ,EOM[i] ) ) : 

_  >  od: 

'  This  section  makes  a  change  of  notation  so  equations  are  compatible  with  standard  MATLAB  notation 

for  state  variables  and  inputs 

[  >  xl :  =  [ ] : xldot :  =  [ ] : 

[>  for  i  from  1  to  vectdim(DOF)  do  xldot:  =  [op(xldot)  ,x[i]  ]  od: 

■>  for  i  from  vectdim  (DOF) +1  to  2*  vectdim  (DOF)  do  xl  :  =  [op  (xl)  ,x  [i]  ] 
od: 

[  >  setX:={ } : 

■  >  for  i  from  1  to  vectdim (DOF)  do 

>  setX : =setX  union  { dDOFq [ i ] =xldot [ i ] } : 

>  setX:=setX  union  {DOFq[i]=xl [i] } : 

_  >  od: 

'  >  setX; 

{<78  =  ^18’  ^^8  =  ^8>  ^6  =  ^\6>  -  ^7’  =  ^17’  ^^9  =  ^9’  %  =  ^19’  ^^10  =  ^10’  ^10= 

=:  Xi5,  dq^  =  rg,  q^  =  dq^  =  X3,  q^  =  dq^  =  x^,  q^  =  dq^  =Xy,qy=  Xjj,  d%  =  x^ } 

■  >  setXl  :={abs  (1  ,x  [1]  )=0  ,abs  (1  ,x[2]  )=0  ,abs  (1  ,x[33  )=0  ,abs  (1  ,x[4  ]  )=0}  ; 

setXl  ■-  { abs(  1,  X3 )  =  0,  abs(  1,  xj  =  0,  abs(  1,  x^)  =  0,  abs(  1,  x^)  =  0 } 

[  >  A1 : =subs (setX , op (A) )  : 

[>  fl  :=subs  (setX,op  (f ) )  : 

[  >  f2 :=subs (setXl , op (f 1) ) : 

[  >  B : =augment (A1 , f 2 ) : 
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’  >  readlib (fortran) ; 

[>  #fortran  (B, optimized) 


proc(x)  ...  end 

This  statement  converts  equations  to  computer  code 
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APPENDIX  B 


Equations  of  Motion  Generated  for  a  Three  Bladed  Simplified  Rotor 

Model  with  MAPLE 


■  The  Mowing  is  an  excerpt  from  the  MAPLE  worksheet  which  was  programmed  to  carry  out  the 

Tjagrangi^n  derivation  of  the  equations  of  motion  for  a  3  bladed  coupled  rotor  fuselage  system. 

- - 

[> 

[  FUSELAGE  X-DIRECTION; 

>  EOM1[1]=0; 

—mb^  cos(%5 )  el  —  mbj  cos(%4)  el  —  wb^  cos(%2)  Q  el  +  A/[  %7 


-  mb^  R  cos(%5)  cos(Ci(0)  *^i(0 

+  2mb^R  sin(%5)  n  sin(Ci(0)  %6  -  mb^  R  cos(%5)  sin(;i(0) 


r^i 


+  m6j  R  sin(%5)  sin((;i(0)  -  2  mb^  R  cos(%5)  H  cos(Ci(0) 


+  mb,  R  sin(%5)  sin(Ci(0)  %6'  +K,  u,{t)-mb,  R  sin(%5)  cosCCjCO) 

cos(%5)  cos(C,(0)  -"*^2-^cos(%4)  cos(C2(0) 

+  2mb^R  sin(%4)  O.  sin(C2(0)  %3  -  mb^  R  cos(%4)  cos(C2(0) 


-  mb^  R  cos(  %4 )  sin(  C2(  0 ) 


V5/ 


^2(0 


,2  ^2 


+  w62i?sin(%4)aSin(C2(0) 


2  mb^  R  cos(%4)  n  cos(<!^2(0)  *^3  +  mb^  R  sin(%4)  sin(^2(0)  °/o3 


-  mb-  R  sin(  %4 )  cos(  C2(  0 ) 


mb.  R  cos(%2)  cos(C3(0) 


+  2mb^R  sm(%2)  n  sitiCCjCO)  %1  -  ^  cos(%2)  cos{C,^{t.))  %1 


-  m/^3  R  cos(%2)  sin(C3(0) 


C 


\di 


2  ^3 


+  mb.  R  siii(%2)  Q  sin((^3(0) 


-  2  mb.  R  cos(%2)  Q  cos{C,^{t))  %1  +  W2a3  i?  sin(%2)  sin(C3(0) 


-m^>3  7?sin(%2)cos(C3(0) 


\dt 


C3(0 


.2  ^3 


=  0 
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%1  :=^C3(0 

%2  :=  /  +  O, 


%3  :--«<) 

%4  O  /  +  $2 
%5  :=  n  r  +  <I)i 


%6:=-C,(0 


%7:-  — »,(() 

L  5r 

[  FUSELAGE  Y-DIRECTION; 

r  >  EOMl [2]=0; 


-ff/63  sin( %6 )  el  —  mb^  sin(  %4 )Q^  el  -  mh^  sin(  %2 )  e/  +  Cj  i  i/2(  0  j  +  ^2  %( 0 


+M2  %7  +  %7  +  w^>2  +  wZ>3  %7  -  mbj_  R  sin(%2)  sin(C2(0)  <^2(0 

-  mb.  R  cos(%2)  sin(C2(0)  -2mb^R  sin(%2)  005(^2(0) 


cos(%2)  sm((;2(0)  %l  +mb^  R  cos(%6)  cos(^;3(0)  [^^3(0 
-  mb.  R  sin(%6)  Cf  cos((;3(0)  -2mb^R  cos(%6)  O  sin(C3(0) 


-  wZ>3  i?  sm(%6)  cos(C3(0)  %5^  -  mb^  R  sin(%6)  sin(C3(0)  [^'^3(0 

-  m&3  R  cos(%6)  sm(C3(0)  -  2  sm(%6)  Q  cos(C,^(t))  %5 

-  mb^  R  cos(%6)  sin(C3(0)  -mb^R  sin(%4)  cosCCjCO) 

-  2  mb,  R  cos(%4)  Q  sin(C,(0)  %3  -  i?  sin(%4)  cos(Ci(0) 


-wiii?sin(%4)sm(Ci(0)  “Ci(0  -  cos(%4)  Q  sin(^i(0) 

\dt  J 

-  2  mb,  R  sin(%4)  Q  cos(Ci(0)  °/o3  -  mb^  R  cos(%4)  sin(Ci(0) 


+  /wZ>,  i?  cos(%4) cos(^i(0)  -  wij i?  sin(%2)  Q  cos(^2(0) 

V.5/  y 


-  2  mb^  R  cos(%2)  Q  sinCCjCO)  cos(%2)  cos(C2(0)  1  ^^2(0 
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“  rnb^  R  sin(%2)  cos((!^2(0)  =  0 

ot 

%2  !“  O  /  + 


%3:--«0 
%4  !—  O  /“  + 
%5  :=-Ut) 

dr^ 
%6  :=  n  /  + 


%7  := 


3/ 


“2(0 


.2  2 


[  ROTOR  BLADE  1  LEAD  LAG: 

>  EOMl [3]=0; 

fd  ) 

Ke,  +  Czeta,  [-Ci(Oj  " 


(  q2 


\dt 


«i(0 


,2  "I 


7 


f 

Va^  y 

Ta"  ^ 

T^'^hiO 

\dt  j 


i?cos(%l)sin(i!;i(0) 
^ 


i?  sm(  %  1 )  cos( Ci(  0 )  + 


f  q2 


-  mb^ 


( 

i?sm(%l)  sin(Ci(0) + 

KOt 


R  cos(%l)  cos(^i(0) 


+  mb^  Cf  el  R  sm(Q^(t))  =  0 

%1  :=  Q  ^ 

[  ROTOR  BLADE  2  LEAD  LAG: 

>  EOMl [4]=0; 

Ke2  ^2(0  +  ^^^^2  ^  sin(C2(0 )  +  ^^^2 

^ 


(&  1 

v3/  ) 

-mb-^ 

i?sin(%l)sin(C2(0) 


+ 


-mb-. 


fa" 

17"^''^ 


R  cos(%l )  cos(C2(0)  - 


f  ^ 


\dt^ 


Mi(0 


R  sin(%l )  cos(;^2(0)  +  Czeta.^ 


\dt‘ 

^  a 


Ml(0 


.2  "1 


«o 


i?cos(%l)  sin(C2(0) 


=  0 


%1  :—  Q  /  +  Oj 

[  ROTOR  BLADE  3  LEAD  LAG: 

[■  >  EOMl  [5]=0; 
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Ke,  C3(0  +  Czeta,  [f  CaCoJ  +  ^'^3  "  "'^3  sm(%l  sinCCjCO) 

/  o2  f  d^  '] 

+  mb.  —th{t)  i?cos(%l)cos(C3(0)-"^^3  7I'^'t(0  i?cos(%l)sin((;3(/)) 

VS/’  J  ^ot  J 

^  s’  1 

-mb.  —,um  ;?sin(%l)cos((;3(/))  +  w63n’e7i?sin(C3(0)  =  0 

V5/’  ; 

%1  !“  ^  +  3^3 
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APPENDIX  C 


OPTIMIZED  CODE  GENERATED  BY  MAPLE  FOR  THE  SIMPLE  THREE 
BEADED  COUPLED  ROTOR-  FUSELAGE  MODEL 


MAPLE  converts  the  elements  of  B,  which  is  an  augmented  matrix  B=  [A  f|,  jfrom  their  symbolic 
representation  into  FORTRAN  code  (or  C  code  if  desired). 

>  fortran (B, optimized) / 
t2  =  rob{l)*R 
t3  =  Omega* t 
t4  =  t3+Phi(l) 
t5  =  cos (t4) 
t6  =  sin(x(8) ) 
t7  =  t5*t6 
t9  =  sin(t4) 
tlO  =  cos (x {8 ) ) 
til  =  t9*tl0 
tl3  =  -t2*t7-t2*tll 
tl4  =  mb{2)*R 
tl5  =  t3+Phi(2) 
tl6  =  cos(tl5) 
tl7  =  sin(x(9) ) 
tl8  =  tl6*tl7 
t20  =  sin(tl5) 
t21  =  cos (x{9) ) 
t22  =  t20*t21 
t24  =  -tl4*tl8-tl4*t22 
t25  =  mb<3)*R 
t26  =  t3+Phi(3) 
t27  =  cos(t26) 
t28  =  sin(x{10) ) 
t29  =  t27*t28 
t31  =  sin(t26) 
t32  =  cos (x (10) ) 
t33  =  t31*t32 
t35  =  -t25*t29-t25*t33 
t42  =  Oraega**2 
t43  =  t42*el 
t45  =  t5*t42 
t48  =  t2*t9 
tSO  =  Omega*t6*x (3) 
t52  =  t5*tl0 
t53  =  x(3)**2 
t56  =  t9*t42 
t59  =  t2*t5 
t61  =  Omega*tlO*x (3) 
t63  =  t9*t6 
t68  =  tl6*t42 

t71  =  -v(l)*x(l)*abs (x(l) )-K(l)*x(6)-c(l)*x(l)+mb{l)*t5*t43+t2*t45 

#*tl0-2*t48*t50+t2*t52*t53-t2*t56*t6+2*t59*t61-t2*t63*t53+inb (2) *tl6 

#*t43+tl4*t68*t21 
t72  =  tl4*t20 
t74  =  Omega*tl7*x (4) 
tie  =  tl6*t21 
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til  =  x(4)**2 

t80  =  t20*t42 

t83  =  tl4*tl6 

t85  =  Omega*t21*x ( 4) 

t87  =  t20*tl7 

t92  =  t27*t42 

t95  =  t25*t31 

t97  =  Omega*t28*x (5) 

t99  =  t27*t32 

tlOO  =  x{5)**2 

tl03  =  t31*t42 

tl06  =  t25*t27 

tl08  =  Omega*t32*x (5) 

tllO  =  t31*t28 

tll3  =  -2*t72*t74+tl4*t76*t77-tl4*t80*tl7+2*t83*t85-tl4*t87*t77+inb 
#(3)*t27*t43+t25*t92*t32-2*t95*t97+t25*t99*tl00-t25*tl03*t28+2*tl06 

#*tl08-t25*tll0*tl00 
tll8  =  -t2*t63+t2*t52 
tl21  =  -tl4*t87+tl4*t76 
tl24  =  t25*t99-t25*tll0 

tl45  =  -c(2)*x{2)-v(2)*x(2)*abs(x(2) )-K(2)*x(7)+t25*t92*t28+2*t95* 
#tl08+mb (1) *t9*t43+t2*t56*tl0+2*t59*t50+t2*tll*t53+t2*t45*t6+2*t48* 
#t61+t2*t7*t53 

tl67  =  mb (2) *t20*t43+t25*t29*tl00+tl4*t80*t21+2*t83*t74+-tl4*t22*t7 
#7+tl4*t68*tl7+2*t72*t85+tl4*tl8*t77+iTib  (3)  *t31*t43+t25*tl03*t32+2*t 
#106*t97+t25*t33*tl00 
tl69  =  R**2 

tl76  =  Ks  (1)  *signura(x  (8) -z) 
tl80  =  Ks  (1)  *signuin{x  (8) +z) 
tl85  =  el*R 
tl89  =  x(8)**2 

tl94  =  u ( 1) -2*Vzeta (l)*x(3)*abs(x(3))-tl76*x(8)/2+tl80*z/2+tl80*x{ 
#8) /2+tl76*z/2-inb (1) *t42*tl85*t6-Ke (1) *x (8) -Kd(l) *tl89*x (8) -Ks (1) *x 
#(8)-Czeta(l)*x(3) 
tl98  =  Ks{2)*signum(x(9)-z) 
t203  =  Ks (2) *signum{x (9) +z) 
t213  =  X ( 9) **2 

t218  =  tl98*z/2-tl98*x(9) /2+t203*z/2+t203*x(9) /2-2*Vzeta (2 ) *x ( 4 ) *a 

#bs{x(4)  )+u(2)-inb(2)*t42*tl85*tl7-Ke(2)*x(9)-Kd(2)*t213*x{9)-Ks(2)* 

#x(9)-Czeta(2)*x(4) 
t220  =  x(10)**2 
t226  =  Ks (3) *signum(x (10) -z) 
t231  =  Ks  (3)  *signuin(x  (10) +z) 

t242  =  -Kd(3) *t220*x (10) -Ke (3) *x ( 10) +u (3) -t226*x (10) /2+t226*z/2+t2 
#31*x ( 10) /2+t231*z/2-2*Vzeta (3) *x (5 ) *abs (x  (5) ) -mb (3) *t42*tl85*t28-C 
#zeta(3)*x(5)-Ks(3)*x(10) 


B(l,l) 

= 

mb ( 1 ) +mb ( 2 ) +mb ( 3 ) +M ( 1 ) 

B(l,2) 

= 

0 

B(l,3) 

tl3 

B(l,4) 

t24 

B(l,5) 

t35 

B(l,6) 

t71+tll3 

B(2,l) 

0 

B(2,2) 

= 

M(2)  +mb  (1)  +rab  (2)  +mb  (3) 

B(2,3) 

=: 

tll8 

B(2,4) 

= 

tl21 

B(2,5) 

tl24 
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B(2,6)  =  tl45+tl67 
B(3,l)  =  tl3 
B(3,2)  =  tll8 
B  (3, 3)  =  mb (1) *tl69 
B(3,4)  =  0 
B(3,5)  =  0 
B(3,6)  =  tl94 
B(4,l)  =  t24 
B(4,2)  =  tl21 
B(4,3)  =  0 
B(4,4)  =  mb(2)*tl69 
B{4,5)  =  0 
B(4,6)  =  t218 
B(5,l)  =  t35 
B(5,2)  =  tl24 
B(5,3)  =  0 
B(5,4)  =  0 
B{5,5)  =  mb(3)*tl69 
B(5,6)  =  t242 


APPENDIX  D 


S-FUNCnON  M-FILE  REPRESENTING  THE  DYNAMICS  OF  THE  SIMPLE 
ROTOR-FUSELAGE  THREE  BEADED  MODEL 


fiinrtiiin  [sys*  xOJ  =  hdo3bA(L*»tt»fl*8>Il»I2jI3»I4>B4^) 

%  function  [sySj  lO]  =  hdo3bA(t^,u4log>Il»I2»J3»I4,K^6) 

% 

%  S-function  arguments: 

%  — - 

%  t  =  time 
%  X  =  state  vector 
%  u  =  input  vector 

%  flag  =  switch  used  by  numerical  integration  (simnlation) 

%  routine  to  access  certain  parts  of  tbe  s-function 

% 

%  S-function  input  parameters: 

%  - - - 

% 

%  II  =  [mb(l)^b(2)»mb(3)^(l)^(2)l 
% 

•/•  12  =  [R,Omega,el^] 

% 

%  13  =  lPbi(l)^bi(2)^bi(3)l 

% 

%  14  =  [c(l),c(2),v(l),v(2), 

%  Caeta(l),Czeto(2),Cieta(3), 

%  Vzeta(l),Vzeta(2),Vzeta(3)] 

% 

%  IS  =  [Ke(l)JKe(2)^e(3), 

%  Kd(l)JKd(2);Kd(3), 

%  Ks(l)4^(2)^3), 

•/.  K(l)^2)l 

% 

%  16  =  [xrXi^Yi^lUr2Ur3i, 

%  ^O^Yi^li^i^i] 

% 

%  S-function  to  represent  dynamics  of  3  bladed  coupled  rotor- 
%  fusdage  modd  wbicb  considers  only  inplane  d^rees  of 
%  freedom,  Le^  x  and  y  translational  fusdage  d^rees  of  freedom 
%  and  lead-lag  rotor  blade  d^rees  of  freedom. 
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% 

%  Explanation  of  variables: 
% - 


-> 

-> 

-> 


% 

%  mb 
•/.M 

%R 
%  el  -> 
%  Omega 
•/.z  -> 

•/.Phi  -5 
•/.c  -> 

%  V  -> 

%Czeta  -> 
%  Vzeta  -> 
-> 


•/•K 

•/•Ke 
•/•Kd 
•/•Ks 
%  xr_ 

•/. 


-> 

-> 

_i-> 
i  -> 


mass  of  blade 
effective  mass  of  fuselage 
distance  from  lead-lag  hinge  to  blade  center  of  mass 
blade  hinge  offset 
->  rotor  speed 
angle  at  which  blade  hits  stops 
blade  phase  angle  w.r.t  azimuth  postion 
fuselage  linear  damping 
fuseh^e  hydraulic  damping 
blade  linear  damping 
blade  hydraulic  damping 

effective  stiffness  of  fusdage  (landing  gear  stiffness) 

blade  dastic  spring  constant 

blade  duffing  spring  constant 

blade  stop  effective  spring  constant 

initial  rate 

initial  displacement 


%  Define  input  parameters 


mb(l)=Il(l);mb(2>=n(2);inb(3)=Il(3);M(l)=Il(4);M(2)=Il(5); 

R=I2(l);Omega=I2(2);el=12(3);z=I2(4); 

Phi(l)=I3(l);Phi(2)=I3(2);Phi(3)=I3(3); 

c(1)=I4<1);c(2)=I4(2);v(1)=I4(3);v(2H4<4); 

CzeU(l)=I4(5);Czeta(2)=I4(6);Cieta(3)=I4(7); 

VzeU(l)=I4<8);Vzeta(2)=I4(9);Vzeta<3)=I4(10); 

Ke(l)=I5(l);Ke(2>=I5(2);Ke(3)=I5(3); 

Kd(l)=15(4);Kd(2)=I5(S);Kd(6>=I5(6); 

K5(1>=I5(7);Ks(2>=I5{8);Ks(3)=I5(9); 

K(1)=I5(10);K(2)=K(11); 

xrXi=I6(l);xrYi=I6(2);xrli=I6(3);ir2i=I6(4);xr3i=I6(5); 

xXf=I6<«);xYi=I6(7);xli=I6(8);x2i=I6(9);x3i=I6(10); 


%  S-fiinction  flag  conditionals 
ifflag  =  0 


sys=[10, 0,10^,0,01; 
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iO=[xrXi^Yi^li,xr2i»xr3i»xXi^Yi^li^i^il; 

dseif  abs(Cbig)  =  1 

%  Formulated  equations  of  motion  optimized  for  minimum  number  of  floating 

%  point  operations. 

t2  =  mb(l)*R; 
t3  =  Omega*t; 
t4  =  t3<-Plii(l); 
tS  =  sin(t4); 
t6  =  cos(x(8)); 
t7  =  t5»t6; 
t9  =  cos(t4); 
tlO  =  sin(x(8)); 
tll=t9ni0; 
tl3  =  -t2*t7-t2*tll; 
tl4»mb(2)*R; 
tl5  =  t3+Phi(2); 
tl6  =  cos(tl^; 
tl7  =  sin(x(9)); 
tl8  =  tl4*tl7; 
t20  =  sin(tl5); 
t21  =  cos(x(9)); 

•  t22  =  t20*t21; 

t24  =  -tl4*tl8-tl4*t22; 
t2S  -  mb(3)*R; 

I26  =  t3+Plu(3); 

t27  =  cos(t26); 

t28  =  sin(x(10)); 

t29  =  t27*t28; 

t31 »  sin(t26); 

t32  =  cos(x(10)); 

t33  =  t31*t32; 

t35  =  -t25*t29-t25*t33; 

t40  =  t2*t9; 

t42  =  Omega*t6*x(3); 

t44  =  t5*tl0; 

t45  =  x(3r2; 

t48  -  Omega'^2; 

t49  =  tl6*t48; 

t52»tl4*t20; 

t54  =  Omega*tl7*x(4); 

t56  =  tl6*t21; 

t57  =  x(4r2; 


95 


t60  =  t20*t48; 

t63  =  tl4*tl6; 

t65  =  Omega*t21*x(4); 

t67  =  t20»tl7; 

m  =  t48*el; 

t73  =  t27*t48; 

t76  =  -V(l)»x(l)»abs(i(l))-c(l)*x(l>f2*t40*t42-t2*t44*t45+tl4»t49*... 

t21-2*t52*t54+tl4*t56*t57-tl4*t60»tl7+2»t63*t65-tl4*t67*t57+mb(3)*... 

t27*t71+t25*t73»t32; 

t77  =  t25»t31; 

t79  =  Omega*t28*x(5); 

t81  =  t27*t32; 

182  =  x(5)^2; 

t85  =  t31*t48; 

t88  =  t25*t27; 

t90  =  Oniega*t32*x(5); 

t92  =  t31*t28; 

t97  =  t9»t48; 

tlOO  =  t2*tS; 

tl02  =  Omega*tlO*x(3); 

tl04  =  I9*t6; 

tl07  =  15*148; 

tll3  =  -2*t77*t79+t25*t81*t82-t25*t85*t28+2*t88*l90-t25*l92*t82+mb(l)^. 
*t9*t71+t2*t97*t6-2*tl00*tl02+t2*tl04*t45-t2*tl07*tlCH-mb(2)*tl6... 
*t71-K(l)*x(6); 
tll8  =  -t2*t44+t2*tl04; 
till  =  -tl4*t67+tl4*t56; 
tl24  =  -t25*t92+t25*t81; 

tl46  =  -v(2)*x(2)*ab8(x(2)>.K(2)*x(7)-c(2)*x(2>fmb(2)*t20*t71+tl4*^ 
t60*t21+2*td3*t54+tl4*t22*t57+iiib(3)*t31*t71+tl4*t49*tl7+2*l52*t65+... 
tl4*tl8*t57+t25*t85*t32; 

tl67  =  2*(88*t79+t25*l33*t82+t25*t73*t28+2*t77*l9(H-t25*t29*t82+mb(l).^ 
•t5*t71+t2*tl07*t6+2*t40*tl02+t2*t7*t45+t2*t97*tl(H-2*tlOO*t42+t2... 
nil*t45; 
tl69  =  R^2; 
tl71  =  x(8)''2; 
tl77  =  i^l)*sign(x(8)>z); 
tl85  =  Ks(l)*sigii(x(8>fz); 
tl89  =  el*R; 

tl96  =  -Kd(l)*tl71*x(8)-Ke(l)*x(8>^-tl77*z/2-tl77*x(8y2-2*Vzeto(l),.. 
*x(3)*abs(x(3)>+tl85*z/2+tl85*x(8y2-mb(l)*t48*tl89*tl(H^u(l)^. 
-Ks(l)*x(8)-CzeU(l)*x(3); 
t203  =  Ks<2)*sigii(x(9)-z); 
t208  =  Ks(2)*sigii(x(9>+z); 
t217  =  x(9)^2; 


% 


t222  =  -2*V2eta(2)*i(4)*abs(x(4))+t203*2/2-t203*x(9)/2+t208*2/2+t208... 
*x(9y2-mb(2)*t48*tl89*tl7-^2)*x(9>-C*eta(2)*x(4)-Ke(2)*x(9)-Kd(2)... 

•t217*x(9)+u(2); 
t224  =  1(10^2; 
t237  =  Ks(3)*sign(x(10)-z); 
t242  =  Ks(3)*sigii(x(10)+z); 

1248  =  -Kd(3)*t224*x(10)-Ke(3)*x(10)-mb(3)»t48*tl89*t28-Ks(3)*x(10)... 
-2*Vxeto(3)*x(5)*abs(x(5))-t237*x(10)/2+t237*z/2+t242*a/2-Czeta(3)... 
*x(5)+u(3)+t242*x(10)/2; 

=  inb(l)+mb(2)+mb(3>+M(l); 

B(U)  =  0; 

B(U)  =  tl3; 

B(M)  =  t24; 

B(l^  =  t35; 

B(l,6)  =  t76+tll3; 

B(2,l)  =  0; 

B(2^)  =  M(2)+mb(lHinb(2)+mb(3); 

B(2^)  =  tll8; 

B(2,4)  =  tl21; 

B(2^  =  tl24; 

B(2,6)  =  tl46+tl67; 

B(3,l)  =  tl3; 

B(3^)=tll8; 

B(3^)  =  mb(l)*tl69; 

B(3,4)  =  0; 

B(3^  =  0; 

B(3,6)«tl96; 

B(4,l)  =  t24; 

B(4^)  =  tl21; 

B(4^)  =  0; 

8(4,4)=^  mb(2)ni69; 

B(4^  =  0; 

B(4,6)  =  t222; 

B(5,l)  =  t35; 

B(5^)  =  tl24; 

B(5^)  =  0; 

B(5^)  =  0; 

B(5^  =  mb(3)*tl69; 

B(5,6)  =  t248; 


%  Calculate  derivatives 
[m,n]=size(B); 


sys=zeros(ly2*m); 

sys(l:5)=Al\n; 

5ys(6;10)=x(l;5); 

%  Output  states 

dseif  abs(fUg)  =  3 

sys=x; 

else 

sys  =  D; 

end 
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APPENDIX  E 


MOVING  BLOCK  ANALYSIS  CODE 

The  following  group  of  MATLAB®  programs  can  be  used  to  perform  a  moving 
block  modal  damping  analysis  on  a  signals  that  are  either  unimodal  or  bimodal.  The 
organization  of  the  code  is  as  follows: 

mbloc  is  the  primary  code  and  calls  maxCm  and  frecur ;  maxf2m  calls  getmax  and  ffl 
(  from  the  MATLAB®  Signal  Processing  Toolbox  function  library)  and  dft .  dampA  is  a 
separate  code  that  is  used  to  curve  fit  the  resulting  moving  block  plot. 

MBLOC: 

function  [logFl,tl4ogF2,t2,omega^N^bl=mbIocpC^r) 

%  function  pogF,t]=niblocpt^r) 

% 

%  MBLOC  calculates  the  magnitude  of  the  discrete 
%  Fourier  transforms  of  Mock  segments  of  a  signal 
%  for  moving  block  damping  analysis.  This  code  is 
%  specificalfy  designed  to  handle  a  signal  widi  1  or 
%  2  dominant  modes. 

% 

%  X  ->  vector  which  contains  the  s^al 
%  sr  ->  sampling  rate  at  which  signal  was  obtained 
%  logFl»2>>  vector  containing  the  natural  logs  of  the 
%  moving  block  function  for  each  successive 

•/•  block 

%  ->  vector  containing  the  times  initializing 

%  each  block 

%  omega  ->  firequency  at  which  the  moving  block  function 
%  is  evaluated 

%  N  ^  signal  leugth  after  it  is  padded  with  zeros 

%  Nb  ->  block  length.  _ 

%  f  ->  frequency  spectrum  ofFFT(Oto  Nyquistfireq) 

%  Pn  ->  power  spectrum  (magnitude  of  FFT) 

%  wl>2  ->  vector  of  frequencies  over  zoomed  frequency 

%  interval 

% 

%  Copyright  (c)  1997  by  Chris  S.  Robinson 
%  All  rights  lesmrved 
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*/•  Call  routine  which  determines  frequency  of  interest 
%  and  block  size  for  evaluation  of  moving  block  function 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


|oniega^yNb4‘»Pzx,wl^bsFl,w2^bsF2]=maxf2m(XySr); 


%  Pad  signal  with  zeros  if  length  is  not  a  power  of  2 
•/•  (this  step  b  done  because  the  fft  routine  contained 
%  in  the  function  maxfum  also  pads  the  original  signal 
%  with  zeros  if  necessary).  If  signal  length  b  a  power 
%  of  2  initially  tiien  thb  step  leaves  the  signal 
%  unaltered. 


l=N-length(X); 

z=zeros(l»l); 

X=IXzl; 


%  Evaluate  the  moving  block  function  along  signal  using  frecnr  and  then  fit 
%  resulting  curve  witti  linear  least  squares  fit 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


[logFl4l]=frecur(X,omega(l)tNb(l)ySr); 

PogF2,t2]=frecur^om^(2)»Nb(2)^r); 

pl=polyrit(tl4ogFl,l); 

p2=po^it(t24ogF24); 

fitl=polyval(pl,tl); 

fit2=polyval(p2yt2); 


%  Plot  the  FFT  resufts,  the  zoomed  DFT  results  and  moving  block  function  versus 
%  time 


subplot(3,14) 
plot(nPxi,V); 
zlabd('frequency  (Hz)'); 
yUbeK'radians^2/Hz'); 
title('Moving  Block  Plot'); 
subplot(3,l>2) 

plot(wlyabsFl/r',w2^bsF2/g'); 
xlabd('frequency  (Hz)'); 
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ylabel('niduuis^2/Hz*); 

sub|rfot(3,l^) 

pIot(tl4ogFl,’r%tl4itl,’^->t24ogF2,*g%t24it2,’g-’); 

xIabdCtime  (sec)*); 

ylabel(’Iog(|F(w)|)*); 

grid 

dampl=-pl(l)/(2*pi*omega(l)) 

damp2=-p2(l)/(2*pi*om^a(2)) 

MAXF2M: 


function  [oniegaJJ,NB,f,Pxx,wl^bsFl,w2,absF2]=niarf2ni(X,sr) 

%  fiinction  omega=inax£2mp^r) 

% 

%  MAXF2N  computes  the  2  dominant  maximum  of  a 
%  bt-modal  signal,  X,  in  the  frequency  domain 
%  by  using  a  fit  for  an  initial  estimate  and  then  refining 
%  the  solution  by  dividing  the  interval  bounded  by 
%  the  nearest  harmonics  to  the  fit  solution  into  subintervals. 

%  The  Fourier  coefficients  are  found  at  each  of  the 
%  frequencies  defined  by  the  subintervals,  and  a  nw 
%  miiTininm  is  found.  The  intervals  nearest  to  the  maximum 
%  are  further  subdivided  and  the  maximum  obtained  b 
%  considered  an  adequate  estimate. 

% 

%  X  •>  Vector  containing  tiie  values  of  the  signal 
%  sr  ->  Mmpling  rate  at  which  the  signal  generated/recorded 
%  om^  ->  the  fkequency  of  tiie  dominant  mode  present  in  the 
%  s^al 

%N  o^iengtii  of  signal  padded  witii  zeros 
%  Nb  ->  length  of  signal  Mode  that  will  be  used  for  the 
%  moving  Mock  analysis 

%Pxx  ->  Power  spectrum  of  s^al  (magnitude  of  FFT) 

%f  ->  frequency  spectrum  of  i<  F1‘ 

%  ->  two  dominant  frequencies  estimated  from  FFT,  rdined 

%  estimate  wiO  be  made  about  tiiese  two  frequencies 

%  absFl,2  ->  Power  spectrum  over  the  zoomed  intervals  about 
%  about  die  estimated  frequencies 

%  omegal,2  ->The  refined  estimates  of  two  dominant  frequencies 
%  NB1,2  ->Block  sizes,  corresponding  to  the  refined  frequency 
%  estimates,  to  be  used  for  the  moving  block 

%  calcualtions 

% 

%  Copyist  (c)  1997  by  Chris  S.  Robinson 
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%  An  rights  reserved 


%  Determine  signal  length  and  the  number  of  points  to  be  added 
%  to  make  that  length  a  power  of  2,  then  take  the  fast  Fourier 
%  transform,  instructing  the  fft  routine  to  pad  the  signal  with 
%  the  proper  number  of  zeros  if  necessary. 


n=length(X); 

N=2^(ceU0og2(n))); 

XF=m(X,N); 

%  Take  the  results  of  the  fft  and  determinethe  power  spectrum 
%  of  the  signal 


m=length(XF); 

Pxr=XF.*conj(XFym; 

nyq=sr/2; 

f=nyq»(0;N/2)/(N/2); 

Pii(N/2+2:N>=D; 

Pn(2:N/2)=2*Pix(2:N/2); 

%  Get  a  first  estimate  of  the  signal  frequencies  by  finding  the 
%  frequencies  corresponding  to  the  spikes  in  the  power 
%  spectrum,  then  take  the  largest  two. 


[maxV4ndV]=getmax(Pxx); 

jmaxVs41=sort(maxV); 

maxVs=fliplr(maxVs); 

i=fliplr(t); 

u=f(indV(i(l))); 

l=f(mdV(K2))); 

dOmega=abs(n-l); 


%  Zoom  in  on  estimated  frequency  and  take  discrete  Fourier 
%  transform  for  frequencies  on  an  interval  around  the  initial 
%  estimates  in  order  to  refine  the  fft  resolution 
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nb»41oor0ength(X)/2); 

kl=nb/sr*u; 

k2’=Db/sr*l; 

h=ceil(length(X)*dOmega/nyq) 

Nb=nb-10*h;nb+10*h; 

wl=Nb.''(-l)*kl*sr, 

w2=Nb.'^(-l)*k2*sr; 

absFl=dft(X^r^,wl); 

absF2*dft^t»sr»Nb»w2)} 


%  Determine  refined  fkvquenqr  estimate  from  results  of  discrete 
%  Fourier  transforms  of  zoomed  interval  by  finding  frequency 
%  corresponds  to  the  maximum  in  the  zoomed  power  spectrum 
%  plot 


[niaxF14ndl]=getmax(absFl); 

[niaxF2»iod2]^ctmax(absF2); 

[dummy4ndla]’=max(maxFl); 

indlb=indl(indla); 

[dummy jnd2a]=max(maxF2); 
ind2b=ind2(ind2a); 


%  Return  the  frequency  and  Mock  size  for  two  modes 

NBl»Nb(indlb); 

om^al=wl(mdlb); 

NB2»Nb(ind2b); 
om^a2^w2{ind2b)y 
NB=[NB1  NB21; 
om^a-[omegal  omega2]; 

FRECUR: 

function  [logF,tl=frecurpi:,om^a^,sr) 

%  function  po^,t]=frecur(X,omega^,sr) 

%  .  . 
%  FRECUR  evaluates  the  moving  Mode  function  for  a  signal 

%  given  the  frequency  of  interest  and  the  Mock  lengtii 

%  using  the  recursion  method. 


103 


%  (only  good  for  boxcar  windowing) 

% 

%  Copyright  (c)  1997  by  Chris  S.  Robinson 
%  All  rights  reserved 


%  Evaluate  the  Fourier  coefficients  for  the  initial  block. 
%  This  step  also  initializes  the  recursion  formula  for 
%  evaluation  of  all  subsequent  blocks. 


N=dength(X); 

kb=omega*Nb/sr; 

Xb=X(l:Nb); 

c=cos(2  *pi*kb/Nb*(0:  Nb-1)); 

s=sin(2*pi*kb/Nb*(0:Nb-l)); 

a(l)=2/Nb*sum(Xb.*c); 

b(l)=2/Nb*sum(Xb.»s); 

t(l>=0; 


%  Evalnte  the  Fourier  coefficients  for  the  reamaining 
%  blocks  by  applying  the  recursion  formula 


forT=2:N-Nb; 

a(l>=(a(T-l)+2/Nb*(X(NbfT-l>-X(T-l)))*cos(2*pi*kb/Nb).~ 

+b(T-l)*sin(2»pi*kb/Nb); 

b(T)=-{a(T-l>f2/Nb*(X(NbfT-l>.X(T-l)))»sin(2*pi*kb/Nb)... 

+b(T-l)*cos(2*pi*kb/Nb); 

t(TKr-l)/sn 

end 


%  Evaluate  the  natural  log  of  the  moving  block  function  for 
%  each  block. 

logF^l/2»log(a.''2+b.^2); 


GETMAX; 


function  [niaxV4ndV]=getmax(X) 
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%  function  [maxV,indV]=getinaxpC) 

% 

%  GETMAX  determines  the  relative  maximum  points 
%  in  a  vector  of  data  (X). 

nH[engdi(X); 
connt=0; 
for  i=l:n-2 
a=X(H-l)-X(i); 
b=X(i+l)-X(i+2); 
if(sign(a)  =  sign(b)  &  a  >  0) 
count=count+l; 
niaxV(count)*=X(i+l); 
mdV(count)^l; 
md 
end 

PFT; 

function  absF^=dft(X^r^,w) 

%  function  absF^ft(X^r^,w) 

% 

%  DFT  computes  the  discrete  Fourier  transform 
%  magnitude  of  a  s^naly  Xy  sampled  at  a  rate*  sr 
%  at  the  frequent  w  (w  is  in  Bx), 

% 

%X  ->  Vector  containing  signal 
%  sr  •>  sampling  rate  at  which  signal  was  created 

%Nb  ^Vector  of  number  of  points  in  sub-Mock  of 
%  signal  over  which  the  discrete  Fourier 
%  transform  will  be  applied 
%  w  ->  Vector  of  frequencies  over  which  the  discrete 
%  Fourier  transform  will  calculated  (each  has  a 
%  corresponding  block  size  fit>m  the  Nb  vector) 

% 

%  Copyright  (c)  1997  by  Chris  S.  Robmson 
%  All  rights  resmved 

%  Evaluate  discrete  Fourier  Transform  by  calculating 
%  the  Fouriw  coefficients  at  the  frequency  of  interest 
%  Return  a  vector  containing  the  mi^nitudes  of  the  dft 
%  at  the  frequencies  contained  in  the  vector  w. 
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for  i=l:iength(w) 


N=Nb(i); 

k=N*w(i)/sr; 

i=0:N-l; 

tl=2*pi*k*i/N; 

c=cos(tl); 

s=siii(tl); 

C=X(l:N).»c; 

S=X(l:N).*s; 

t=2/N*suiii(C); 

b=2/N*suiii(S); 

*bsF(i)=a"^2+b^2; 


end 

PAMPA: 


function  3PKljunpA(logF,t,tstart,tstop,oniega) 

*/•  function  y=danipA(logF,t)tstart,tstop, omega) 

•/. 

%  function  DAMPA  perfomu  least  squares  fh  of  moving 
%  block  data  and  calcuahes  the  damping  modal  damping 
%  firom  tile  slope  of  the  fh. 

•/. 

%  logF  ->  vector  of  natural  logs  of  the  moving  block 
%  function  values. 

%  t  ->  time  vector  which  corresponds  to  the  times 
%  at  the  beginning  of  the  each  block  in  the 

%  moving  Mock  analysis. 

*/•  tstart  ->  Start  time  for  section  of  moving  block  plot 
%  to  be  least  squares  fitted. 

%  tstop  ->  Stop  time  for  section  of  moving  block  plot 
%  to  be  least  squares  fitted. 

%  omega  ->  frequency  of  used  in  moving  block  analysis. 
% 

%  Copyright  (c)  1997  by  Chris  S.  Robinson 
%  AO  rights  reserved 

•/. 


%  Extract  the  pertinent  section  of  the  moving  block  plot 
%  according  to  the  user  defined  start  and  stop  times 


u=ones(14cngth(t)); 

utl=t5tart*a; 

ut2=tstop*u; 

[dummyyindl]=min(abs(t-utl)); 

[dummyyind2]=min(abs(t-ut2)); 

%  Perform  a  first  order  polynomial  fit  to  moving  block  plot 


p^lyfit(t(indl:ind2)4ogF(indl:ind2),l); 


%  Use  the  resulting  slope  firom  the  least  squares  fit  to 
%  detmmine  the  modal  damping  and  return  this  value 


y=-p(l)/(2*pi*omega); 
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APPENDIX  F 


input  file  for  complex  ROTOR-FUSELAGE  MODEL 

%  This  m-file  serves  as  input  file  for  running  the  simulink 
%  S-function  hdo3B.m. 

% 

•/••/••/••/•%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/* 
%  Hdo  Physical  and  Aerodynamic  parameters 
%%%%%%%%%%%%%%%%%%%%%%%*/.%%%%%%%%*/»%%% 

%  Distance  from  fiisdage  center  of  mass  to  hub  length). 

h=.7907; 

%  Hinge  offset  length). 
el=^791; 

%  Length  of  rotor  blade  (length). 

R»  2.3809; 


%  Mass  of  rotor  blades  (mass). 

mb(l)=0.01432;  %  LBS/g  (SLUGS) 

mb(2)=0.01432; 

mb(3)=0.01432; 

%  Effective  mass  of  fusdage  and  moments  of  inertia  (mass 
%  or  mass*length^2). 


Ml^l.5486; 

M2=1.30d0; 

Ill=.13497; 

I22=.4669; 

112=0; 

%  Blade  azimuth  phase  angles  (radians) 

Phi(l)=0; 

Plii(2)==2*pi/3; 

Phi(3)=4»pi/3; 

%  Rotor  speed  (radians  per/sec) 
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Omeg*=75^9; 

%  Blade  spring  stiffnesses 

%  Linear  springs  for  lead-lag  (moment/radian) 

kll(l)=143.85; 

kll(2>=143.S5; 

Ul(3)=143.85; 

%  Linear  springs  for  flap  (moment/radian) 

kn(2)=31J8; 

kn(3)=3138; 

%  Doffing  springs  for  lead-lag  (moment/radian^3) 

kl3(l)=«; 

kl3(2>=0; 

kl3(3)=0; 

%  Duffing  springs  for  flap  (moment/radian^S) 

ld3(l)=0; 

ld3<2>=0; 

kl3<3)=0; 

%  Blade  damping  constants 

%  Damping  in  lead-lag  (moment/(rad/scc)) 

cl(l>=3^; 

cl(2)=3^; 

d(3>=3^; 

%  Damping  in  flap  (moment/(rad/sec)) 

cf(l)=0; 

cf(2)=0; 

cf(3)=0; 

%  Fusdage  effective  stiffness 


no 


%  Translational  (force/length) 


KT1=3000; 

KT2=3000; 

%  Rotational  (moment/radian) 

KR1=28035; 

KR2=46.60; 

%  Fusdage  damping  constants 

%  Translational  linear  (force/(longth/sec)) 

CT1=0; 

CT2=0; 

%  Translational  nonlinear  (force/(lengtli/sec)^2) 

VT1=0; 

VT2=0; 

%  Rotational  linear  (moment/(rad/sec)) 

CR1=1.061449; 

CR2=1^9852; 

%  Rotational  nonlinear  (moment/(rad/sec)^2) 

VR1=0; 

VR2=0; 

%  Aerodynamic  parameters. 

%  lift  cnrve  slope  (1/radian) 
a=2*pi; 

%  Parasite  drag  coefficimit 
cd0=0.0079; 

%  Air  density  (mass/length^S) 
riiol=0.002377; 
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%  Rotor  chord  (length) 
c=0.137; 

%  Advance  ratio 
niu=0; 

%  Inflow  ratio 

lambda=0;  %  set  to  zero  for  aU  cases  until  trim 
%  routine  is  setup 

%  Initial  conditions 

%  Fuselage  translational  rates  (Icngth/sec) 

xrtXi=0; 

xrtYi=0; 

%  Fuselage  rotational  rates  (radians/sec) 

irrXi=.l; 

xitYfO; 

%  Blade  lead-lag  rates  (radians/sec) 

xrllM); 

ir21i=0; 

ir31i=0; 

%  Blade  flap  rates  (radians/sec) 

xrlfi=0; 

ir2fi=0; 

xr3fi=0; 

%  Fuselage  translational  displacements  (length) 

xOUfC; 

xtYi=0; 

%  Fuselage  rotational  displacements  (radians) 
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irXi=0; 

irYi=0; 

%  Blade  lead-lag  displacement  (radians) 

xlli=0; 

x21^; 

x31i=0; 

%  Blade  flap  displacement  (radians) 

xlfr*©; 

x2fi=0; 

x3fi=0; 

%  Form  input  matrices  from  above  parameters 

Ils[byelyR*mb(l))ml>(2))mb(3)yM13^*If 

Omega^lii(l)^lii(2)4’lii(3)]; 

I2^[kll(l)^l(2)Jdl(3)4d3(l)4d3(2)4d3(3);.. 

kfl(l)4rfl(2)4rfl(3)Jd3(l)4rf3(2)JsO(3);.^ 

cf(l),cf(2),cl(3),cl(l),cl(2),cl(3)J; 

13=[KTl,KT2,KR14at2;.. 

Cn,CT2,CRl,CR2;.^ 

VT1,VT2,VR1,VR21; 

I4=(ayCd0^ol,Cyma4ambda] ; 

I5=[xrtXi^Yi^Xi^YU**lli^2fi^3fi4rlli^21i^31il; 

I<»*[xtXi^Yi^Xi^Yi^lfi^fi,x3fi^lU^li^b1; 
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